<jupyter_start><jupyter_text>Initialization<jupyter_code>import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import glob
import time
from datetime import datetime, timedelta
# Plotting initialization
import seaborn as sns
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd

%matplotlib inline

# Initialize style
plt.style.use("acm_ieee_latex_pubstyle.txt")

SINGLE_COLUMN = SINGLE_WIDTH = 8.0
DOUBLE_COLUMN = DOUBLE_WIDTH = 16.0

def width_height(width=SINGLE_COLUMN, height=None, columns=1, rows=1):
    """Width driven sizing for seaborn multi-grid layouts
    
    Derive height & aspect kwargs for facetgrid, pairgrid & friends with a target sizing
    in terms of width and height.  Convenient for creating figures that has
    to be fit into a certain width, while maintaining the correct aspect ratios of the fonts.
    Call this function to derive the height & aspect
    
    Currently, ACM & IEEE latex, single column figure defaults to a width of 8inch
    and the double column figure defaults to a width of 8 inch for the correct text scaling
    
    Usage:
    - sns.FacetGrid( ... , **width_height())
    - sns.PairGrid( ..., **width_height(columns=X, rows=Y))

    :param width: target width in inches (4 inch for a single column figure is default)
    :param height: target heigth in inches (width * 3/4 - golden ratio is the default)
    :param columns: target # of columns for facetgrid
    :param rows: target # of rows for facetgrid
    """
    if height is None:
        height = width * 3/4
    ratio = float(width) / (float(height) / float(rows))
    return {"height": (float(width) / ratio), "aspect": ratio / float(columns)}
# Dask setup block
import os
import pwd
import glob
import pandas as pd
from distributed import LocalCluster, Client
import dask
import dask.dataframe as dd

# Cleanup
try:
    client.shutdown()
    client.close()
except Exception as e:
    pass

# Configuration
LOCALDIR = "/tmp/dask"
dask.config.set({'worker.memory': {'target': False, 'spill': False, 'pause': 0.8, 'terminate': 0.95}})

# Cluster creation
cluster = LocalCluster(processes=True, n_workers=2, threads_per_worker=15,
                       dashboard_address=":8787", local_directory=LOCALDIR,
                       memory_limit="16GB")

client = Client(cluster)
cluster
client<jupyter_output>/opt/conda/lib/python3.8/site-packages/distributed/node.py:151: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 46779 instead
  warnings.warn(<jupyter_text>Dataset Selection and Test Read<jupyter_code># Component temperature & CEP data
TEMP_PATH = "/gpfs/alpine/stf218/proj-shared/data/lake/summit_thermal_cluster_comptype_cep"
# Power data
POWER_PATH = "/gpfs/alpine/stf218/proj-shared/data/lake/summit_power_temp_openbmc/power_ts_10s_2"

# Valid dates for temperature data
# - 20200724 ~ 20201124
sdate = pd.to_datetime("2020-07-24")
edate = pd.to_datetime("2020-09-30")
#edate = pd.to_datetime("2020-08-30")
date_keys_with_temp = [dt.date().strftime("%Y%m%d") for dt in list(pd.date_range(sdate, edate - timedelta(days=1),freq='d'))]

# Work to do
date_keys_with_temp = [
    work for work in date_keys_with_temp
        if os.access(f"{TEMP_PATH}/{work}.csv", os.F_OK) and os.access(f"{POWER_PATH}/{work}.csv", os.F_OK)
]

# Power only data
# - 20200101 ~ 20201124
sdate = pd.to_datetime("2020-01-01")
edate = pd.to_datetime("2020-12-31")
date_keys_only_power = [dt.date().strftime("%Y%m%d") for dt in list(pd.date_range(sdate, edate - timedelta(days=1),freq='d'))]

# Work to do
date_keys_only_power = [
    work for work in date_keys_only_power
        if os.access(f"{TEMP_PATH}/{work}.csv", os.F_OK) and os.access(f"{POWER_PATH}/{work}.csv", os.F_OK)
]
TIMESTAMP = "timestamp"
# Power consumption
VALID_COUNT = "count_inp"
NODE_COUNT = "size_inp"
POWER = "sum_inp"
MAX_POWER = "max_inp"
MIN_POWER = "min_inp"
STD_POWER = "std_inp"

# Component temperature
GPU_CORE_MEAN = "gpu_core.mean" 
GPU_CORE_MAX = "gpu_core.max"
GPU_CORE_MIN = "gpu_core.min"
CPU_CORE_MEAN = "cpu_core.mean"
CPU_CORE_MAX = "cpu_core.max"
CPU_CORE_MIN = "cpu_core.min"

# Secondary loop
SUPPLY_TEMP = "mtwst"
RETURN_TEMP = "mtwrt"
SUPPLY_FLOW = "mtwflw"

# Efficiency
PUE = "k100_pue"
KW_PER_TON = "cep_kw_per_ton"
MTW_TON = "cep_mtw_tons"
CT_TON = "cep_cooling_tower_tons"
CHW_TON = "cep_chilled_water_tons"
WETBULB_TEMP = "cep_outside_air_wet_bulb_temp"

# Primary loop flow
CT_FLOWRATE = "ct_water_flowrate"
CHW_FLOWRATE = "chw_flowrate"
MAKEUP_FLOW = "cep_make_up_flow"

# List
CLUSTER_POWER = [VALID_COUNT, NODE_COUNT, POWER, MAX_POWER, MIN_POWER, STD_POWER]
COMPONENTS = [GPU_CORE_MEAN, GPU_CORE_MAX, CPU_CORE_MEAN, CPU_CORE_MAX]
SECONDARY_LOOP = [SUPPLY_TEMP, RETURN_TEMP, SUPPLY_FLOW]
PRIMARY_LOOP = [CT_FLOWRATE, CHW_FLOWRATE, MAKEUP_FLOW]
EFFICIENCY = [PUE, KW_PER_TON, MTW_TON, CT_TON, CHW_TON]
WEATHER = [WETBULB_TEMP,]<jupyter_output><empty_output><jupyter_text>Test Data  Power data<jupyter_code>df = pd.read_csv(POWER_PATH + f"/{date_keys_with_temp[0]}.csv")
df.head()
df.columns<jupyter_output><empty_output><jupyter_text>Temperature data<jupyter_code>df = pd.read_csv(TEMP_PATH + f"/{date_keys_with_temp[0]}.csv")
df.head()
df.columns<jupyter_output><empty_output><jupyter_text>Detect movements and carve out segments from both power & tempPrepare a large table that has many short segments* Select only the dates that have temperature data as we lack temperature sensor data for certain dates* For each dates, join power & temp data, sort  - Find rows where the power delta is larger than the specified criteria and store the indexes  - For each indexes, lookup the original time series, carve out the rows around the center row      - When carving out:          - filter out starting points based on the specified criteria - we should bin them with 1MW intervals          - filter out segments that has too many missing nodes          - translate the indicies into offsets that starts from 0      - The timestamp of the center row should be in a category      - Also, categorize based on the specified criteria<jupyter_code>%%time

import functools
import operator
    
# Number of timesteps to take before and after
SEGMENT_FRONT =12
SEGMENT_REAR = 64


PERIODS = 4


def handle_day(date_key, direction="rise"):
    """Process the day"""
    BAND_WATTS = (1000 * 1000)
    
    thresholds = [
        {"class": 13, "start": 13000000, "end": 14000000},
        {"class": 12, "start": 12000000, "end": 13000000},
        {"class": 11, "start": 11000000, "end": 12000000},
        {"class": 10, "start": 10000000, "end": 11000000},
        {"class": 9, "start": 9000000, "end": 10000000},
        {"class": 8, "start": 8000000, "end": 9000000},
        {"class": 7, "start": 7000000, "end": 8000000},
        {"class": 6, "start": 6000000, "end": 7000000},
        {"class": 5, "start": 5000000, "end": 6000000},
        {"class": 4, "start": 4000000, "end": 5000000},
        {"class": 3, "start": 3000000, "end": 4000000},
        {"class": 2, "start": 2000000, "end": 3000000},
        {"class": 1, "start": 1000000, "end": 2000000},
    ]

    # Filenames
    temp_file = f"{TEMP_PATH}/{date_key}.csv"
    power_file = f"{POWER_PATH}/{date_key}.csv"

    # Join file
    power_df = pd.read_csv(power_file, usecols=["timestamp"]+ CLUSTER_POWER).sort_index()
    temp_df = pd.read_csv(temp_file).sort_index()
    all_df = pd.merge(power_df, temp_df, on="timestamp")
    max_index = all_df['timestamp'].count()
    
    # Get diff
    diff_df = all_df[POWER].diff(periods=PERIODS)

    # Inner carving function
    def carve_edge(row, threshold):
        """Segment carving"""
        # Carve out the segment
        index = int(row["index"])
        diff = row[POWER]
        
        # Filter incomplete segments
        if index < SEGMENT_FRONT or index >= max_index - SEGMENT_REAR:
            return []

        segment = all_df.iloc[index - SEGMENT_FRONT:index + SEGMENT_REAR].reset_index(drop=True)
        start_inp_class_mw = int(segment.iloc[SEGMENT_FRONT][POWER] / BAND_WATTS)

        # Filter segments with insufficient data
        ts = segment[segment[NODE_COUNT] < (4626 - 200)]['timestamp'].count()
        if ts > 0:
            return []
        
        # Categorize the segment
        segment['edge_index'] = index
        segment['segment_timestamp'] = segment['timestamp']
        segment['start_class_mw'] = start_inp_class_mw
        segment['amplitude'] = diff
        segment['amplitude_class_mw'] = threshold['class']
        return segment

    # Iterate through the thrshold spec
    def carve_threshold(threshold):
        """Carve out a threshold based segment set"""
        # setup mask based on rise / fall
        assert(direction in ['rise', 'fall'])
        mask = ((diff_df >= threshold['start']) & (diff_df < threshold['end']))
        if direction != "rise":
            mask = ((- diff_df) >= threshold['start']) & ((- diff_df) < threshold['end'])

        # Apply mask and filter out errors
        edge_indexes = diff_df[mask].reset_index()
            
        # Apply segment carving and filter out empty batches
        segment_lst = [
            seg for seg in list(edge_indexes.apply(carve_edge, args=(threshold,), axis=1))
                if type(seg) is pd.DataFrame
        ]
        return segment_lst
    
    # Iterate through thresholds
    segment_list = [carve_threshold(thr_spec) for thr_spec in thresholds]
    segment_list = [seg for seg in segment_list if len(seg) > 0]
    if not segment_list:
        return []

    # Done with the segments, concatenate them and exit
    flat_list = functools.reduce(operator.iconcat, segment_list, [])
    return flat_list

def filter_duplicates(flat, direction="rise"):
    """Filter out duplicates"""
    segments = sorted(flat, key=lambda x: -x.iloc[0]['amplitude_class_mw'])
    indices = []
    ret_flat = []
    
    def compare(x):
        start = index - SEGMENT_FRONT
        end = index + SEGMENT_REAR
        res = (x >= start) and (x < end)
        return res
    
    # Iterate segments and accumulate indices
    for i, segment in enumerate(segments):
        index = segment.iloc[0]['edge_index']
        existing_indicies = list(filter(compare, indices))
        if not existing_indicies:
            ret_flat.append(segment)
            indices.append(index)
    return ret_flat

# Distribute
def get_segments(date_keys, direction="rise"):
    """Get segments from the full dataset"""
    import dask.bag as db
    segments = [seg for seg in db.from_sequence(date_keys).map(handle_day, direction=direction).compute() if seg]
    
    item_count = [len(seg) for seg in segments]
    print(f"Segments collected: {direction} {sum(item_count)}, {item_count}")
    
    flat = functools.reduce(operator.iconcat, segments, [])
    filtered = filter_duplicates(flat, direction=direction)
    print(f"Segments survived: {direction} {len(filtered)}")
    return pd.concat(filtered)


def get_segments_slow(date_keys, direction="rise"):
    segments = [handle_day(date_key, direction=direction) for date_key in date_keys]
    
    item_count = [len(seg) for seg in segments]
    print(f"Segments collected: {direction} {sum(item_count)}, {item_count}")
    
    flat = functools.reduce(operator.iconcat, segments, [])
    filtered = filter_duplicates(flat, direction=direction)
    print(f"Segments survived: {direction} {len(filtered)}")
    return pd.concat(filtered)

# Test process a particular day
segment_list = handle_day("20200728")
%%time
# Get all the rise & falling segments
rise_df = get_segments(date_keys_with_temp, direction="rise")
fall_df = get_segments(date_keys_with_temp, direction="fall")

rise_df['type'] = "rise"
fall_df['type'] = "fall"
rf_df = pd.concat([rise_df, fall_df])
rf_df<jupyter_output>Segments collected: rise 942, [22, 260, 9, 31, 18, 10, 2, 3, 25, 8, 13, 29, 50, 24, 10, 33, 27, 76, 3, 2, 7, 12, 25, 6, 28, 2, 1, 2, 8, 24, 13, 12, 5, 4, 14, 1, 7, 104, 3, 5, 4]
Segments survived: rise 156
Segments collected: fall 958, [23, 245, 8, 4, 29, 31, 7, 24, 4, 28, 2, 2, 73, 10, 49, 24, 1, 68, 7, 2, 13, 18, 41, 10, 21, 5, 1, 1, 12, 27, 15, 12, 9, 1, 4, 1, 1, 2, 105, 6, 6, 6]
Segments survived: fall 152
CPU times: user 3.07 s, sys: 224 ms, total: 3.3 s
Wall time: 17.2 s<jupyter_text>Power Dynamics - Summer  Summer power dynamicsThis one will be the substitute for the rise / fall edge power graphs.  Should give a more comprehensive view.Will be using a 16 x 4 strip for this one<jupyter_code>%%time
# Get all the rise & falling segments
rise_df = get_segments(date_keys_with_temp, direction="rise")
fall_df = get_segments(date_keys_with_temp, direction="fall")

rise_df['type'] = "rise"
fall_df['type'] = "fall"
rf_df_summer = pd.concat([rise_df, fall_df])
%%time
# initialize
rf_df = rf_df_summer
plt.close()

# Override acm-ieee-latex-pubstyle
plt.figure(figsize=(16, 2))
sns.despine(top=False, right=False)

# Layout
gs = mpl.gridspec.GridSpec(2, 7, height_ratios=[1, 1]) 
gs.update(hspace=0.1, wspace=0.1)

plt.xlabel('Time (minutes)', loc='center')

##################################################
# Comparison
def plot(rf_df, ax, amp_class, direction, metric, kwargs):
    mask = (
        (rf_df['type'] == direction) & (rf_df['amplitude_class_mw'] == amp_class)
    )
    df = rf_df[mask].reset_index()
    # TODO: Some massage (melting)
    sns.lineplot(ax=ax, data=df, x="index", y=metric,
                 sort=True, ci=95, **kwargs)
    
layout = [
    # Rise
    (0, gs[0, 0], 1, '1MW', 'rise', POWER),
    (1, gs[0, 1], 2, '2MW', 'rise', POWER),
    (2, gs[0, 2], 3, '3MW', 'rise', POWER),
    (3, gs[0, 3], 4, '4MW', 'rise', POWER),
    (4, gs[0, 4], 5, '5MW', 'rise', POWER),
    (5, gs[0, 5], 6, '6MW', 'rise', POWER),
    (6, gs[0, 6], 7, '7MW', 'rise', POWER),
    # Fall
    (7, gs[1, 0], 1, '1MW', 'rise', PUE),
    (8, gs[1, 1], 2, '2MW', 'rise', PUE),
    (9, gs[1, 2], 3, '3MW', 'rise', PUE),
    (10, gs[1, 3], 4, '4MW', 'rise', PUE),
    (11, gs[1, 4], 5, '5MW', 'rise', PUE),
    (12, gs[1, 5], 6, '6MW', 'rise', PUE),
    (13, gs[1, 6], 7, '7MW', 'rise', PUE),
]
for i, loc, amp_class, amp_label, direction, metric in layout:
    # Subplot item
    ax = plt.subplot(loc)
    kwargs = {
        "linewidth": 1,
        "legend": False,
    }
    # Count number of segments per amplitude class
    seg_count = rf_df[rf_df['type'] == 'rise'].groupby('amplitude_class_mw')['timestamp'].count()

    # Location customization
    if i not in [0, 7]:
        plt.setp(ax.get_yticklabels(), visible=False)
    if i in [0, 1, 2, 3, 4, 5, 6]:
        ax.set_title(f"{amp_label} - {seg_count[amp_class]:,}")
        plt.setp(ax.get_xticklabels(), visible=False)

    # Plot
    plot(rf_df, ax, amp_class, direction, metric, kwargs)
    
    # X ticks
    xtck = [0, 12, 24, 36, 48]  # ticks
    xlab = [-1, 0, 1, 2, 3]   # minutes
    ax.set_xlim([0, SEGMENT_REAR])
    ax.set_xticks(xtck)
    ax.set_xticklabels(xlab)
    ax.set_xlabel("")
    
    if i in [0, 1, 2, 3, 4, 5, 6]:
        # Y ticks - Power in MW
        MW = 1000000
        ytck = [4 * MW, 6 * MW, 8 * MW, 10 * MW]
        ytlb = [4, 6, 8, 10]
        ylim = [2000000, 12000000]
        ylab = "Power\n(MW)"
    else:
        # Y ticks - PUE
        ytck = [1.0, 1.1, 1.2, 1.3, 1.4]
        ytlb = [1.0, 1.1, 1.2, 1.3, 1.4]
        ylim = [1.0, 1.4]
        ylab = "PUE"
        
    ax.set_ylim(ylim)
    ax.set_yticks(ytck)
    ax.set_yticklabels(ytlb)
    ax.set_ylabel(ylab)

    # Overlay
    ax.axvline(SEGMENT_FRONT, color='r', linewidth=0.3)

fig = plt.gcf()
fig.text(0.5, -0.09, 'Time after Rising Edges (mins)', ha='center')
    
plt.subplots_adjust(hspace=.0)
plt.savefig("../plots/power_dynamics_per_amp_summer.pdf", format="pdf", bbox_inches='tight', pad_inches=0.01)
plt.show()<jupyter_output><empty_output><jupyter_text>Power Dynamics - Full Year  Full year power dynamicsThis one will be the substitute for the rise / fall edge power graphs.  Should give a more comprehensive view.Will be using a 16 x 4 strip for this one<jupyter_code>%%time
# Get all the rise & falling segments
rise_df = get_segments(date_keys_only_power, direction="rise")
fall_df = get_segments(date_keys_only_power, direction="fall")

rise_df['type'] = "rise"
fall_df['type'] = "fall"
rf_df_fullyr = pd.concat([rise_df, fall_df])
%%time
# initialize
rf_df = rf_df_fullyr
plt.close()

# Override acm-ieee-latex-pubstyle
plt.figure(figsize=(16, 2))
sns.despine(top=False, right=False)

# Layout
gs = mpl.gridspec.GridSpec(2, 7, height_ratios=[1, 1]) 
gs.update(hspace=0.1, wspace=0.1)

plt.xlabel('Time (minutes)', loc='center')

##################################################
# Comparison
def plot(rf_df, ax, amp_class, direction, metric, kwargs):
    mask = (
        (rf_df['type'] == direction) & (rf_df['amplitude_class_mw'] == amp_class)
    )
    df = rf_df[mask].reset_index()
    # TODO: Some massage (melting)
    sns.lineplot(ax=ax, data=df, x="index", y=metric,
                 sort=True, ci=95, **kwargs)
    
layout = [
    # Rise
    (0, gs[0, 0], 1, '1MW', 'rise', POWER),
    (1, gs[0, 1], 2, '2MW', 'rise', POWER),
    (2, gs[0, 2], 3, '3MW', 'rise', POWER),
    (3, gs[0, 3], 4, '4MW', 'rise', POWER),
    (4, gs[0, 4], 5, '5MW', 'rise', POWER),
    (5, gs[0, 5], 6, '6MW', 'rise', POWER),
    (6, gs[0, 6], 7, '7MW', 'rise', POWER),
    # Fall
    (7, gs[1, 0], 1, '1MW', 'rise', PUE),
    (8, gs[1, 1], 2, '2MW', 'rise', PUE),
    (9, gs[1, 2], 3, '3MW', 'rise', PUE),
    (10, gs[1, 3], 4, '4MW', 'rise', PUE),
    (11, gs[1, 4], 5, '5MW', 'rise', PUE),
    (12, gs[1, 5], 6, '6MW', 'rise', PUE),
    (13, gs[1, 6], 7, '7MW', 'rise', PUE),
]
for i, loc, amp_class, amp_label, direction, metric in layout:
    # Subplot item
    ax = plt.subplot(loc)
    kwargs = {
        "linewidth": 1,
        "legend": False,
    }
    # Count number of segments per amplitude class
    seg_count = rf_df[rf_df['type'] == 'rise'].groupby('amplitude_class_mw')['timestamp'].count()

    # Location customization
    if i not in [0, 7]:
        plt.setp(ax.get_yticklabels(), visible=False)
    if i in [0, 1, 2, 3, 4, 5, 6]:
        ax.set_title(f"{amp_label} - {seg_count[amp_class]:,}")
        plt.setp(ax.get_xticklabels(), visible=False)

    # Plot
    plot(rf_df, ax, amp_class, direction, metric, kwargs)
    
    # X ticks
    xtck = [0, 12, 24, 36, 48]  # ticks
    xlab = [-1, 0, 1, 2, 3]   # minutes
    ax.set_xlim([0, SEGMENT_REAR])
    ax.set_xticks(xtck)
    ax.set_xticklabels(xlab)
    ax.set_xlabel("")
    
    if i in [0, 1, 2, 3, 4, 5, 6]:
        # Y ticks - Power in MW
        MW = 1000000
        ytck = [4 * MW, 6 * MW, 8 * MW, 10 * MW]
        ytlb = [4, 6, 8, 10]
        ylim = [2000000, 12000000]
        ylab = "Power\n(MW)"
    else:
        # Y ticks - PUE
        ytck = [1.0, 1.1, 1.2, 1.3, 1.4]
        ytlb = [1.0, 1.1, 1.2, 1.3, 1.4]
        ylim = [1.0, 1.4]
        ylab = "PUE"
        
    ax.set_ylim(ylim)
    ax.set_yticks(ytck)
    ax.set_yticklabels(ytlb)
    ax.set_ylabel(ylab)

    # Overlay
    ax.axvline(SEGMENT_FRONT, color='r', linewidth=0.3)

fig = plt.gcf()
fig.text(0.5, -0.09, 'Time after Rising Edges (mins)', ha='center')

plt.subplots_adjust(hspace=.0)
plt.savefig("../plots/power_dynamics_per_amp_fullyr.pdf", format="pdf", bbox_inches='tight', pad_inches=0.01)
plt.show()<jupyter_output><empty_output>