<jupyter_start><jupyter_text>Table of ContentsUsing Principal Component Analysis (PCA) on Deschutes Fermentation DataThe GoalIssues to considerThis NotebookPart 1: Define System ParametersPart 2: Load &amp; Clean DataPart 3: Define Utility FunctionsPart 4: Process and Summarize DataPart 5: Apply PCA to FeatureDFPart 6 (BONUS): Apply clustering to Principal Components:Your TaskPart 1. Define System Parameters1a. Define relevant main variables1b. Define auxiliary variables1c. PCA coordinates selectionPart 2. Load &amp; Clean Data2a. Initialize Hub client2b. Download data from OCS with Data ViewsFind predefined multi-asset Data View for Fermenters 31-36 and PCAData View StructureGet Interpolated Data from Data Views2c. Clean all_brands_dfPart 3. Define Utility FunctionsPart 4. Process and Summarize DataPart 5. Apply PCA to FeatureDF5a. Define and scale the feature space5b. Apply PCA to the feature space5c. Find number of principal components sufficent to explain data5d. Visualize Deschutes data using the principal components5e. What's in an outlier?Part 6. Apply clustering to Principal Components (BONUS)6a. Apply k-means clustering on the principal components6b. Plot principal components and color by kmeans cluster --- Using Principal Component Analysis (PCA) on Deschutes Fermentation Data---It is often useful to be able to quickly determine whether and which of the production batches are out-of-spec. The ability to determine the origin or the cause of a process failure or a bad batch is valuable to industry.In this notebook, we will use Principal Component Analysis (PCA) to determine which batches are out-of-spec without having to analyze the time series data for every attribute of erach beer batch, of which there are dozens. We will then use the properties of PCA to identify what may have caused a batch to be out-of-spec.First, make sure the latest version library for OCS is installed.<jupyter_code>try:
    import ocs_academic_hub
except ModuleNotFoundError:
    !pip install ocs-academic-hub==1.0.1
    import ocs_academic_hub

!pip install plotly_express==0.4.1 sklearn<jupyter_output><empty_output><jupyter_text>Then import all the required packages<jupyter_code># Academic Hub module 
from ocs_academic_hub.datahub import hub_login

# For OCS configuration and data manipulation
import configparser
from dateutil import parser
import datetime as dt

# analysis packages
import pandas as pd
import numpy as np

# packages for PCA. sklearn is standard Python library for machine learning
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans

# packages for plotting
from plotly.subplots import make_subplots
import plotly.graph_objs as go
print("--- import completed ---")<jupyter_output><empty_output><jupyter_text>--- The GoalUse Principal Component Analysis (PCA) to identify the bad beer batches in 300 days worth of process data involving 4 beer brands (Realtime Hops, Trois Lacs, Grey Horse and LIT), brewed across 6 fermentor vessels (31-36). Also, use the properties of Principal Components (PCs or "scores") to identify what features are most responsible for contributing to the out-of-spec qualities of outlier beer batches. Issues to consider* Real process data is dirty (sensors can fail)* The stages may have been mislabeled (e.g. entries in "Cooling" phase may have been mislabled as being in "Diacetyl Rest")* Different cooling temperatures for each beer batch, and even for each zone in a fermentor vessel for a given beer batch* There can be various volumes of beer in the fermentors* Data can have abrupt, aphysical jumps in temperature --- This Notebook Part 1: Define System ParametersSpecify parameters of interest like brands, fermentors, time period to analyze, time granularity, attributes. Part 2: Load & Clean Datause OSIsoft Cloud Services (OCS) to obtain process data from Deschutes Brewery. Part 3: Define Utility FunctionsCreate utility functions for calculating the time elapsed per batch, identifying fermentation stages without relying on the "Status" labels, removing the bad batches, and summarizing the data. Part 4: Process and Summarize DataExecute the functions created in Part 3 for each brand in each fermentor vessel. Collect the summary of the data into the dataframe `FeatureDF`. Part 5: Apply PCA to `FeatureDF``FeatureDF` is a high-dimensional dataset that is not easily human-processable. Use Principal Components Analysis (PCA) to obtain new components composed of linear combinations of the most "relevant" input features in `FeatureDF`( [5a](section_5a)-[5c](section_5c)); use those components to graphically visualize `FeatureDF` and identify the bad batches [5d](section_5d); then identify the input features which contributed to the out-of-spec characteristics of the outlier batches ([5e](section_5e)). Part 6 (BONUS): Apply clustering to Principal Components:Use the k-means clustering algorithm to classify data. --- Your TaskSome cells have missing pieces of code and are marked with `TODO` items in comments. Complete all the missing code to obtain a working notebook. If the notebook runs as inteded, using the default seetings in [Part 1a](section_1a), one should obtain this plot in [Part 5e](section_5e):![](https://academichub.blob.core.windows.net/images/nb3_pca_scatterplot_new_v6.png)(NOTE: `TODO` items can be found in [Part 3](section_3) only, but this exercise will be easier if you go through it in sequential order). --- Part 1. Define System Parameters---  1a. Define relevant main variablesWe want 3 years of brewing data from 6 fermentor vessels (31-36), beginning on January 11, 2017.<jupyter_code># fermenters to analyze as a list of integers
FERMENTORS_OF_INTEREST = [str(i) for i in range(31, 37)]

# brands of interest (brand `None` indicates that no beer in actually processing)
BRANDS_OF_INTEREST = ["Realtime Hops", "Trois Lacs", "Grey Horse", "LIT"]


# relevant days over which to analyze fermentors - 3 years
TRAINING_DAYS = 3 * 365

# time period in which to analyze fermentation data
START_INDEX = "2017-01-11T00:00"
END_INDEX = (parser.parse(START_INDEX) + dt.timedelta(days=TRAINING_DAYS)).isoformat()

# define time grandularity for data. in this case, get events at 15 minute intervals
INDEX_INTERVAL = "00:15:00"
print("--- main variables defined ---")<jupyter_output><empty_output><jupyter_text>1b. Define auxiliary variablesMake code more readable by avoiding to type long strings repeatedly.<jupyter_code># TIC PV column names (temperatures)
TIC_PV_COLUMNS = ["Bottom TIC PV", "Middle TIC PV", "Top TIC PV"]

# TIC OUT column names (cooling intensities)
TIC_OUT_COLUMNS = ["Bottom TIC OUT", "Middle TIC OUT", "Top TIC OUT"]

# list of relevant attributes. these attributes should not have any bad inputs
ATTRIBUTES_OF_INTEREST = (
    ["Brand", "Status", "Volume", "ADF"] + TIC_PV_COLUMNS + TIC_OUT_COLUMNS
)


# define stages involved in overall fermentation process
ALL_BREWING_STAGES = [
    "Filling",
    "Fermentation",
    "Free Rise",
    "Diacetyl Rest",
    "Cooling",
    "Maturation",
    "Ready to Transfer",
    "Emptying",
    "Clean",
    "CIP",
]

# define stages which can actually affect beer quality
IMPORTANT_BREWING_STAGES = ALL_BREWING_STAGES[
    1:6
]  # Fermentation->Free Rise->Diacetyl Rest->Cooling


# Stages to consider as part of the fermentation stage
FERMENTATION_STAGES = ["Fermentation"]  #  "Free Rise", "Diacetyl Rest"]

# assign an integer value for each brand
brand_dictionary = {brand: (i + 1) for i, brand in enumerate(BRANDS_OF_INTEREST)}
print("--- auxiliary variables defined ---")<jupyter_output><empty_output><jupyter_text>1c. PCA coordinates selection<jupyter_code># ============================================================
# PCA COORDINATES SELECTION

# For temperature sensors with multiple positions, the list
# position_active selects which ones are selected for PCA.
#
# NOTE: at least one position should be True otherwise no sensor
# will be used in PCA even if a 1 is present for the first item in
# PCList below
#
#                  TOP,   MIDDLE, BOTTOM
position_active = [False, False, True]

# List of the names of coordinates to be used for PCA.
# The first item of each indicates if the coordinates is used for PCA (1 for active).
# For sensors with position, see position_active above.
#
PClist = [
    (1, ["Brand"]),  # brand identifier - FIRST COLUMN SHOULD ALWAYS BE 1
    (1, ["Batch"]),  # batch identifier - FIRST COLUMN SHOULD ALWAYS BE 1
    (1, ["FERM_ADF"]),  # variance of ADF during primary fermentation
    (0, ["COOL_ADF"]),  # average adf after primary fermentation
    (0, ["Volume"]),  # average volume after primary fermentation
    (1, ["FERM_Duration"]),  # total time PRIMARY FERMENTATION lasted
    # average value for cooling intensity during PRIMARY FERMENTATION at different zones
    (1, ["FERM_avg_coolint_top", "FERM_avg_coolint_middle", "FERM_avg_coolint_bottom"]),
    # average value for temperature during PRIMARY FERMENTATION at different zones
    (1, ["FERM_avg_temp_top", "FERM_avg_temp_middle", "FERM_avg_temp_bottom"]),
    # maximum value for temperature during PRIMARY FERMENTATION at different zones
    (1, ["FERM_max_temp_top", "FERM_max_temp_middle", "FERM_max_temp_bottom"]),
    # minimum value for temperature during PRIMARY FERMENTATION at different zones
    (1, ["FERM_min_temp_top", "FERM_min_temp_middle", "FERM_min_temp_bottom"]),
    # variance of the temperature during PRIMARY FERMENTATION at different zones
    (1, ["FERM_var_temp_top", "FERM_var_temp_middle", "FERM_var_temp_bottom"]),
    # average of the variance between zones per event during PRIMARY FERMENTATION
    (1, ["FERM_var_between_zones"]),
    (0, ["COOL_Duration"]),  # total time COOLING lasted
    # average value for cooling intensity during COOLING at different zones
    (0, ["COOL_avg_coolint_top", "COOL_avg_coolint_middle", "COOL_avg_coolint_bottom"]),
    # average value for temperature during COOLING at different zones
    (0, ["COOL_avg_temp_top", "COOL_avg_temp_middle", "COOL_avg_temp_bottom"]),
    # maximum value for temperature during COOLING at different zones
    (0, ["COOL_max_temp_top", "COOL_max_temp_middle", "COOL_max_temp_bottom"]),
    # minimum value for temperature during COOLING at different zones
    (0, ["COOL_min_temp_top", "COOL_min_temp_middle", "COOL_min_temp_bottom"]),
    # initial value for temperature during COOLING at different zones
    (0, ["COOL_init_temp_top", "COOL_init_temp_middle", "COOL_init_temp_bottom"]),
    # variance of the temperature during PRIMARY FERMENTATION at different zones
    (0, ["COOL_var_temp_top", "COOL_var_temp_middle", "COOL_var_temp_bottom"]),
    # average of the variance between zones per event during COOLING
    (0, ["COOL_var_between_zones"]),
]


# ==================== END PCA COORDINATES SELECTION ================
print("--- PCA coordinates selection done ---")<jupyter_output><empty_output><jupyter_text>--- Part 2. Load & Clean Data---The Deschutes Brewery data needs to be retrieved from OSIsoft Cloud Services (OCS). We need to retrieve data for the `FERMENTORS_OF_INTEREST` from time period `START_INDEX` to `END_INDEX` at a resolution of `INDEX_INTERVAL` between each event. The data retrieved from OCS is loaded into a dataframe, which is then filtered to ensure it will not contain a missing value for any value in the `ATTRIBUTES_OF_INTEREST` before it is processed.  2a. Initialize Hub client**Execute the cell below and follow the indicated steps to log in (an AVEVA banner would show up)**<jupyter_code>widget, hub_client = hub_login()
widget
dataset = hub_client.current_dataset()
namespace_id = hub_client.namespace_of(dataset)
print(f"Current dataset is {dataset} in namespace_id: {namespace_id}")<jupyter_output><empty_output><jupyter_text>2b. Download data from OCS with Data Views**Development tip (WARNING: executing dataview takes up to 3 minutes)** Development of a notebook involves running code over and over, so you'll want to avoid long running steps when possible. This is why you can run the cell below once, with the resulting dataframe saved in variable `all_brands_df`. If you don't change any of its input parameters, `all_brands_df` is still valid and can be reused.   Find predefined multi-asset Data View for Fermenters 31-36 and PCA<jupyter_code>dv_id = hub_client.asset_dataviews(filter="pca", asset="FV31", multiple_asset=True)[0]
print(dv_id)<jupyter_output><empty_output><jupyter_text>Data View Structure<jupyter_code>print(hub_client.dataview_definition(namespace_id, dv_id).to_string(index=False))<jupyter_output><empty_output><jupyter_text>Get Interpolated Data from Data Views<jupyter_code># retrieve data from OCS
# TODO: complete call to multiple data views
# =========== STUDENT BEGIN ==========

all_brands_df = @@@ Your code here @@@.@@@ Your code here @@@_pd(...,...,...,...,INDEX_INTERVAL)

# =========== STUDENT END ==========

# preview `all_brands_df`
all_brands_df.head()
# dataframe summary
all_brands_df.info()<jupyter_output><empty_output><jupyter_text>2c. Clean `all_brands_df`Drops all entries which have missing value in any columns in listed in `ATTRIBUTES_OF_INTEREST`.<jupyter_code># drop rows with NaNs in attributes of interest
all_brands_df = all_brands_df.dropna(subset=ATTRIBUTES_OF_INTEREST)

all_brands_df["Asset_Id"] = all_brands_df["Asset_Id"].apply(lambda vessel: vessel[-2:])

# all_brands_df.describe(include='all')
all_brands_df.info()
# PCA will be performed on the 4 brands with the most fermentations
# Note 1: those are upper bounds since a fermentation can be discarded if incomplete or with bad data
all_brands_df.groupby(by=["Fermentation Start Time"]).first()["Brand"].value_counts()<jupyter_output><empty_output><jupyter_text>--- Part 3. Define Utility Functions---Given a dataframe containing the fermentation data for one brand, first find all the different fermentation batches; then for each batch, we need to identify which stage each event belongs to.However, the "Status" labels are unreliable, so the stages have to be indetified through other means. Since it can be difficult to identify the transition between different stages, we consider only two rough categories of stages: (1) Pre (not) cooling and (2) Cooling.The  general "Precooling" stage will contain the literal "Fermentation", "Free Rise", "Diacetyl Rest" stages, while the general "Cooling" stage will contain the literal "Cooling" and "Maturation" stages. To distinguish "Precooling" and "Cooling" rows, a new boolean column "Cooling_Stage" contains a False or True value respectively.   We then want to create a summary of each batch, such that we can capture the general property of each batch. The summary for each batch will be collected in the dataframe `FeatureDF`.| Function | Description ||----------|-------------||`fermentation_times`| Identify when fermentation begins; offset timestamps to begin on fermentation starts; add unique label to each fermentation batch||`remove_batches_without_phase` | Certain stages are critical to flavor profile of beer; batches without these stages are just bad||`organize_by_cooling`| finds when cooling starts and then relabels statuses as being either in 'Precooling' or 'Cooling' stages. Necessary because entries may have been mislabelled.||`summarize_by_batch`| summarizes the time evolution of 10 attributes into 41 scalars. the summary is stored in a dataframe `ProcessedData`|The utility functions can all be found in the cell below (don't forget to execute):<jupyter_code>class NoCoolingStagesFound(Exception):
    pass


def fermentation_times(brand_df, fermentation_starts, shift=0):
    """
    DESCRIPTION: Identify when fermentation begins and &
        offsets timestamps to start when fermentation begins.
        Also adds a label to each new fermentation batch
    RETURNS: the input dataframe with offset timestamps, labels for each batch
    """

    # create new df columns "tsf" and "batch". NOTE: "tsf" --> "time since fermentation"
    # TODO: initialize "tsf" with 100000 and "batch" with -1
    # =========== STUDENT BEGIN ==========
    brand_df = brand_df.assign(**{"tsf": 100000, @@@ Your code here @@@})
    # =========== STUDENT END ==========

    # reset index for fermentation starts, so that position index for .iterrows() is correct
    fermentation_starts.reset_index(inplace=True)

    # calculate the time elapsed for each fermentation step belonging to a new fermentation stage
    for count, ferm_begin in fermentation_starts.iterrows():

        # get the timestamp for each start of fermentation
        fermentation_start_time = pd.Timestamp(ferm_begin["Timestamp"])

        # get booleans for new fermentation steps, where a new fermentation step is any step which exceeds the 
        # new fermentation time
        new_fermentation_steps = (
            brand_df["Timestamp"].apply(lambda t: pd.Timestamp(t))
            >= fermentation_start_time
        )

        # give each new fermentation step a new integer label
        brand_df.loc[new_fermentation_steps, "batch"] = count + shift

        # find time elapsed since start of fermentation, in units of float days
        brand_df.loc[new_fermentation_steps, "tsf"] = brand_df.loc[
            new_fermentation_steps, "Timestamp"
        ].apply(
            lambda t: (pd.Timestamp(t) - fermentation_start_time).total_seconds()
            / 86400.0
        )

    # include only steps with non-negative "batch" labels and "tsf" not exceeding 100000 days
    offset_df = brand_df[(brand_df["tsf"] <= 100000) & (brand_df["batch"] >= 0)]

    # find the new label shift
    shift = offset_df["batch"].max()

    return offset_df, shift
print("--- fermentation_times defined ---")
def no_batch_without_phase(offset_df, ref_df):

    # find all batches represented in `ref_df`; these are valid batches
    valid_batches = ref_df["batch"].unique()

    # we only want batches in `valid_batches`
    # TODO: Complete the filter expression
    # =========== STUDENT BEGIN ==========
    offset_df = offset_df[@@@ Your code here @@@@.isin(@@@ Your code here @@@)]
    # =========== STUDENT END ==========

    return offset_df


def organize_by_cooling(offset_df, fermentation_starts):
    """
    DESCRIPTION: finds when cooling starts and then relabels statuses as being either in
        'Precooling' or 'Cooling' stages. Necessary because entries may have been mislabelled.
    RETURNS: relabeled offset_df
    """

    # find events which are stictly in cooling stage. When cooling intensity (TIC OUT) exceeds 99.99 in
    # the three vessel zones, the batch is in the cooling stage.
    # TODO: Complete the expression
    # =========== STUDENT BEGIN ==========
    boolean = (
        (offset_df["Top TIC OUT"] > 99.99)
        & (@@@ Your code here @@@)
        & (@@@ Your code here @@@)
    )
    # =========== STUDENT END ==========

    # filter for events which are strictly for cooling
    strict_cooling_df = offset_df.loc[boolean, :]

    # the start of cooling is the first instance of a strict cooling event for a batch (given by 'batch')
    cooling_starts = strict_cooling_df.groupby("batch").first().reset_index()

    # filter out batches with no cooling stages so it won't be used for PCA
    offset_df = no_batch_without_phase(offset_df, cooling_starts)

    # if there are no batches with cooling stages, don't continue with data processing
    if offset_df.empty:
        raise NoCoolingStagesFound
        return None, None
    else:
        print(f"\t{len(cooling_starts)} Cooling starts")

    # relabel batch statuses as either 'pre-cooling' or 'cooling'
    for __, row in cooling_starts.iterrows():

        # label for current batch
        current_batch = row["batch"]

        # when cooling starts for current batch
        cooling_start_time = row["tsf"]

        # get booleans for current batch
        boolean = offset_df["batch"] == current_batch

        # Add a new "Cooling_Stage" column which is False or True for "Precooling" or "Cooling" respectively
        # TODO: Complete the expression
        # HINT: "Cooling_Stage " should be False if tsf < cooling start time, True otherwise
        # =========== STUDENT BEGIN ==========
         offset_df.loc[boolean, "Cooling_Stage"] = offset_df.loc[boolean, "tsf"].apply(
            lambda t: @@@ Your code here @@@ if @@@ Your code here @@@  else True
        )
        # =========== STUDENT END ==========

    return offset_df
print("--- organize_by_cooling defined ---")
# Goal of the following definitions is from PCList, which is a list of pairs of the following format:
#    (boolean, list of principal components)
# to generate a list of booleans called PCA_df_selected which length is the concatenation of all
# list of principal components in PCList. If a list in PCList has length > 1, the boolean value is
# is copied for all element in the list
#
flatten = lambda l: [item for sublist in l for item in sublist]
PClist_unzip = [list(t) for t in zip(*PClist)]
PC_columns_name = flatten(PClist_unzip[1])
PCA_df_selected = []
[PCA_df_selected.extend([x] * len(y)) for x, y in PClist]


def summarize_by_batch(offset_df):
    """
    DESCRIPTION: summarizes the time evolution of 10 attributes into 41 scalars.
    RETURNS: dataframe containing summary of processed data for one brand in a
        fermentor vessel.
    """

    # get all unique batch labels available in offset_df
    all_batches = offset_df["batch"].unique()

    # get number of batches available in offset_df
    num_batches = len(all_batches)

    # initialize a pandas dataframe with all elements as NaN (Not a Number)
    PCA_df = pd.DataFrame(
        data=np.zeros((num_batches, len(PC_columns_name))), columns=PC_columns_name
    )
    PCA_df[:] = np.nan

    # fill PCA_df with summary of time series data per batch
    for i, batch in enumerate(all_batches):

        # divide offset_df for some given batch into cooling vs primary fermentation stages
        # ferm_df has rows belonging to fermentation
        #
        ferm_df = offset_df[
            (offset_df["batch"] == batch)
            & (offset_df["Status"].isin(FERMENTATION_STAGES))
        ]

        # cool_df has rows belonging to cooling (i.e. "Cooling_Stage" is True) 
        # and to the current `batch` for ferm_df above
        #
        # TODO: Complete the filter expression
        # =========== STUDENT BEGIN ==========
        cool_df = offset_df[
            (@@@ Your code here @@@ == batch) & (@@@ Your code here @@@ == True)
        ]
        # =========== STUDENT END ==========

        # store batch, ADF, and column info in PCA_df
        # select on active components according to PCList
        #
        PCA_df.iloc[i, 0] = brand_dictionary[ferm_df["Brand"].iloc[0]]
        PCA_df.iloc[i, 1] = batch
        if PCA_df_selected[2]:
            PCA_df.iloc[i, 2] = ferm_df["ADF"].mean()  # average ADF in precooling stage
        if PCA_df_selected[3]:
            PCA_df.iloc[i, 3] = cool_df["ADF"].mean()  # average ADF in cooling stage
        if PCA_df_selected[4]:
            # average Volume in cooling stage
            PCA_df.iloc[i, 4] = cool_df["Volume"].mean()

        if PCA_df_selected[5]:
            # summaraize data in primary fermentation stages
            PCA_df.iloc[i, 5] = ferm_df["tsf"].max()  # duration of precooling stages
        if PCA_df_selected[6]:
            # average values of temperature and cooling intensity for each zone during precooling
            PCA_df.iloc[i, 6:9] = [
                ferm_df[c].mean(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_OUT_COLUMNS)
            ]
        if PCA_df_selected[9]:
            # average values of temperature and cooling intensity for each zone during precooling
            PCA_df.iloc[i, 9:12] = [
                ferm_df[c].mean(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]
        if PCA_df_selected[12]:
            # maximum temperature for each zone during precooling
            PCA_df.iloc[i, 12:15] = [
                ferm_df[c].max(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]
        if PCA_df_selected[15]:
            # minimum temperature for each zone during precooling
            PCA_df.iloc[i, 15:18] = [
                ferm_df[c].min(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]
        if PCA_df_selected[18]:
            # variance of the temperature of each zone cureing precooling
            PCA_df.iloc[i, 18:21] = [
                ferm_df[c].var(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]
        if PCA_df_selected[21]:
            # average of the variance of the temperatures in each zone during precooling
            PCA_df.iloc[i, 21] = (ferm_df[TIC_PV_COLUMNS].var(axis=1)).mean()

        # summarize data in cooling stages
        if PCA_df_selected[22]:
            # duration of cooling stages
            PCA_df.iloc[i, 22] = cool_df["tsf"].max() - cool_df["tsf"].min()
        if PCA_df_selected[23]:
            # average values of temperature and cooling intensity for each zone during cooling
            PCA_df.iloc[i, 23:26] = [
                cool_df[c].mean(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_OUT_COLUMNS)
            ]
        if PCA_df_selected[26]:
            # average values of temperature and cooling intensity for each zone during cooling
            PCA_df.iloc[i, 26:29] = [
                cool_df[c].mean(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]
        if PCA_df_selected[29]:
            # maximum temperature for each zone during cooling
            PCA_df.iloc[i, 29:32] = [  # list(cool_df[TIC_PV_COLUMNS].max(axis=0
                cool_df[c].max(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]
        if PCA_df_selected[32]:
            # minimum temperature for each zone during cooling
            PCA_df.iloc[i, 32:35] = [  # list(cool_df[TIC_PV_COLUMNS].min(axis=0))
                cool_df[c].min(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]
        if PCA_df_selected[35]:
            # variance of the temperature of each zone during cooling
            PCA_df.iloc[i, 35:38] = [  # list(cool_df[TIC_PV_COLUMNS].var(axis=0))
                cool_df[c].var(axis=0) if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]
        if PCA_df_selected[38]:
            # average of the variance of the temperatures in each zone during cooling
            PCA_df.iloc[i, 38] = (cool_df[TIC_PV_COLUMNS].var(axis=1)).mean()
        if PCA_df_selected[39]:
            # initial temperature in cooling stage
            PCA_df.iloc[i, 39:42] = [  # list(cool_df[TIC_PV_COLUMNS].iloc[0])
                cool_df[c].iloc[0] if position_active[ci] else np.nan
                for ci, c in enumerate(TIC_PV_COLUMNS)
            ]

    return PCA_df
print("--- summarize_by_batch defined ---")<jupyter_output><empty_output><jupyter_text>--- Part 4. Process and Summarize Data ---Run the cell below to call the previously defined utility functions and  create the dataframe `FeatureDF` which contains a summary of the time series data from `all_brands_df`. The columns of `FeatureDF` can be used as input features for Principal Component Analysis (PCA).<jupyter_code># initialize dataframe which will be later saved into .csv format to be analyzed by PCA
FeatureDF = pd.DataFrame(columns=PC_columns_name)

# initialize batch labels
label_shift = 1

generator = ((f, b) for f in FERMENTORS_OF_INTEREST for b in BRANDS_OF_INTEREST)
for fermentor, brand in generator:
    print(f"\n===> FV {fermentor}. Brand: {brand}")

    # create reduced df containing only data from relevant vessel and brand
    brand_df = all_brands_df[
        (all_brands_df["Asset_Id"] == fermentor) 
        & (all_brands_df["Brand"] == brand)
    ]
    
    # identify fermentation starts
    fermentation_starts = brand_df[
        (brand_df["Status"] == "Fermentation")
        & (brand_df["Status"].shift(1) != "Fermentation")
    ]
    
    # dont' consider data not in important brewing stages
    brand_df = brand_df[brand_df["Status"].isin(IMPORTANT_BREWING_STAGES)]
    
    # if fermentation did not occur, there is nothing to analyze
    N_fermentations = len(fermentation_starts)
    if N_fermentations == 0:
        print ("\tNo Fermentation data to analyze!")
        continue
    else:
        print (f"\t{N_fermentations} Fermentation starts")
        
    # offset timestamps for brand_df, so t=0 when fermentation starts. 
    offset_df, shift = fermentation_times(brand_df, fermentation_starts, shift=label_shift)
    
    # also assign new labels for each brew batch
    label_shift = shift+1
    
    # remove batches with no maturation phases
    offset_df = no_batch_without_phase(offset_df, offset_df[offset_df["Status"] == "Maturation"])
    
    # identify cooling starts; relabel relabel all Statuses before "Cooling" as "Pre-cooling"
    # if cooling stages do not exist, pca not needed to determine whether batch is bad
    try:
        offset_df = organize_by_cooling(offset_df, fermentation_starts)
    except NoCoolingStagesFound:
        print ("\tNo Cooling data to analyze!")
        continue
        
    # summarize time series data into compact dataframe
    FeatureDF = FeatureDF.append(summarize_by_batch(offset_df))
    
# clean indices and remove unused columns 
FeatureDF = FeatureDF.reset_index()
FeatureDF = FeatureDF.drop(columns="index")
FeatureDF = FeatureDF.dropna(axis=1)
        
# write `FeatureDF` to .csv data file
FeatureDF.to_csv("FeatureDF.csv", sep=",", index=False)

# preview processed data
print("\n\nFinished data processing. `FeatureDF` looks like:")
FeatureDF.head(5)<jupyter_output><empty_output><jupyter_text>--- Part 5. Apply PCA to `FeatureDF`---  5a. Define and scale the feature spaceThe feature space `X` contains all the input features (i.e. columns) in `FeatureDF` which can potentially affect the flavor profile of the beers. When PCA is applied on `X`, we obtain new set of coordinates ("principal components") which are orthogonal to each other and point in the directions of maximum variance (eigenvectors).To prevent any one feature dominating the variance calculations, it is recommended that the feature space `X` be scaled. In this case, the feature space `X` is scaled by "centering" each feature and normalizing them by their standard deviation. Using sklearn, this can be accomplished using `StandardScaler()`.<jupyter_code># define the feature space (variables which affect beer quality)
X = FeatureDF.iloc[:, 2:]

# M = number of observations (rows in X),  N = number of components in the feature space (columns in X)
[M, N] = np.shape(X)

# save unscaled feature space (will only be use for plotting purposes)
Xold = X

# scale the feature space
X = StandardScaler().fit_transform(X)
print("--- feature space defined ---")<jupyter_output><empty_output><jupyter_text>**>>>>>===================THIS NEXT IMMEDIATE CELL CAN BE IGNORED===================**Run the cell below to graphically compare the difference between scaled and unscaled features.<jupyter_code># plotting function
def create_box_plot_object(feature_matrix, N_features, title, legend):
    # print("cols=", feature_matrix.columns)
    plot_objects = []
    for i in range(N_features):
        plot_obj = go.Box(y=np.array(feature_matrix)[:, i], name=legend[i])
        plot_objects.append(plot_obj)

    fig = go.Figure(data=plot_objects, layout=go.Layout(title=title))
    fig.show()
    return

# create a plotly object for the unscaled feature space
create_box_plot_object(Xold, N, "Unscaled Features", list(Xold.columns))

# create a plotly object for the scaled feature space
create_box_plot_object(X, N, "Scaled Features", list(Xold.columns))<jupyter_output><empty_output><jupyter_text>**===========================END IGNORABLE CELL===========================<<<<<**  5b. Apply PCA to the feature spaceWhen PCA is applied to the old feature space `X`, we can obtain a covariance (or correlation) matrix and subsequently eigenvectors and eigenvalues.The eigenvectors, or "Principal Axes", point in the directions of the new feature space, `T`, and in the directions of maximum variance of the old feature space `X`.The eigenvalues,`Lambda`, correspond the the magnitude of the Principal Axes. The eigenvalues reflect the amount of variance explained by the Principal Axes.The new feature space `T` is a matrix of "scores" or "Principal Components". Each Principal Component is a linear combination of the input features in `X`.The coefficients for these linear combinations are the elements of their corresponding Principal Axes or eigenvectors.<jupyter_code># obtain maximum number of principal components
n_comps = min([M, N])

# initialize PCA and get covariance matrix
covariance_matrix = PCA(n_components=n_comps)
covariance_matrix.fit(X)

# get matrix of principal components (a.k.a PCs or scores)
T = covariance_matrix.transform(X)

# get eigenvectors
P = covariance_matrix.components_

# get eigenvalues
Lambda = covariance_matrix.explained_variance_
print("--- cell ok ---")<jupyter_output><empty_output><jupyter_text>**>>>>>===================THIS NEXT IMMEDIATE CELL CAN BE IGNORED===================**Visualize the first $n$ components and how they are related to the input features (columns in `FeatureDF`). A principal component is a linear combination of the input features, and the coefficients are the elements of its eigenvector.<jupyter_code># appropriate number of princiapl components to use
first_n_components = 3

# initialize plots
fig = make_subplots(
    rows=first_n_components, cols=1, shared_xaxes=True, vertical_spacing=0.03
)

for i_component in range(first_n_components):
    plot_object = go.Bar(
        x=FeatureDF.columns.tolist()[2:],
        y=P[i_component, :],
        name=f"PC{i_component+1}",
    )

    fig.append_trace(plot_object, i_component + 1, 1)
    fig["layout"][f"yaxis{i_component+1}"].update(title="Coefficients")
fig["layout"].update(width=800, height=800, xaxis=dict(tickangle=45))
fig.show()<jupyter_output><empty_output><jupyter_text>**===========================END IGNORABLE CELL===========================<<<<<**  5c. Find number of principal components sufficent to explain dataThe eigenvalues of the covariance matrix is equivalent to the variance of a component. We want $l$ components which explain most of the variance in the data. A simple and widely used criterion is the cumulative percent variance (CPV).$ CPV(l) = 100 \left( \frac{\sum_{j=1}^l \lambda_j}{\sum_{j=1}^M \lambda_j} \right) \%$<jupyter_code># compute the variance explained by each component, normalize it, and take the cumulative sum
cumulative_variance = covariance_matrix.explained_variance_ratio_.cumsum()

# plotly object for cumulative variance data
data = go.Scatter(
    x=[item + 1 for item in range(len(cumulative_variance))],
    y=[100.0 * item for item in list(cumulative_variance)],
    mode="lines+markers",
)

# plotly object for reference CPV
lines = []
for var in [75, 87, 99]:
    lines.append(
        go.Scatter(
            x=[0] + [item + 1 for item in range(len(cumulative_variance))],
            y=[var] * (len(cumulative_variance) + 1),
            mode="lines",
            line=dict(color="red", dash="dash"),
        )
    )

# plot layout
layout = go.Layout(
    xaxis=dict(title="Number of Principal Components"),
    yaxis=dict(title="Cumulative Variance (%)"),
    height=500,
    showlegend=False,
)

# create and show plot
fig = go.Figure(data=[data] + lines, layout=layout)
fig.show()
from IPython.display import display, Markdown

display(Markdown("**Question: How many PCs do you think you should use?**"))
# =========== STUDENT BEGIN ==========
display(Markdown(
"""
@@@ Your answer here @@@
"""))
# =========== STUDENT END ==========<jupyter_output><empty_output><jupyter_text>5d. Visualize Deschutes data using the principal componentsRegardless of the appropriate number of principal components from CPV, we are limited to visualizing at most three principal components.<jupyter_code># initialize list for plot objects
obj_list = []

# create plot objects for each brand
for i_brand, brand in enumerate(BRANDS_OF_INTEREST):

    # get which entries in ProcessedData correspond to rband
    boolean = FeatureDF["Brand"] == brand_dictionary[brand]

    # create plotly object
    plot_obj = go.Scatter3d(
        x=T[boolean, 0],
        y=T[boolean, 1],
        z=T[boolean, 2],
        text=list(FeatureDF.loc[boolean, "Batch"]),
        hoverinfo="text",
        hovertemplate="Batch %{text:2.0}",
        name=brand + f" ({boolean.sum()})",
        mode="markers",
        marker=dict(size=7, opacity=0.75),
    )

    # append to list of plotly objects
    obj_list.append(plot_obj)
print("--- cell ok ---")
# create layout
layout = go.Layout(
    scene=dict(
        xaxis=dict(title="PC 1"), yaxis=dict(title="PC 2"), zaxis=dict(title="PC 3")
    ),
    margin=dict(l=0, r=0, b=0, t=0),
    legend=dict(x=1.0, y=0.5),
    height=600,
)

# plot!
fig = go.Figure(data=obj_list, layout=layout)
fig.update_layout(
    title_text=f"<br>PCA scatterplot for the 4 brands {BRANDS_OF_INTEREST}<br> with most data",
    title_x=0.5,
)
fig.show()<jupyter_output><empty_output><jupyter_text>5e. What's in an outlier?We can use **contribution plots** to identify which input features are responsible for making an outlier an outlier. One can create a contribution plot for some given PC and some given observation (i.e. a beer batch). We want to compare how the contributions from an outlier differ from the contributions of in-spec batches. To do this, we need the batch label(s) for the outlier(s) and in-spec/reference batches.A contribution of an input feature to a PC for a given observation is the value of the input feature at that observation times its coefficient (element of eignevector for the PC).To find the batch labels, hover over the data points in the 3-D scatter plot found in Part 5d.<jupyter_code># pick a brand to analyze
brand = "Grey Horse"

# identify the batch label for outlier point(s) - has to be a list!
outlier_batches = [184]

############### >>>>> DO NOT CHANGE SECTION - START >>>>> #############
# filter out irrelevant brands
all_brand_batches = FeatureDF[FeatureDF["Brand"] == brand_dictionary[brand]]

# identify labels for in-spec batches
inspec_batches = [
    x for x in all_brand_batches["Batch"].unique() if x not in outlier_batches
]


def extract_feature_values(ref_lbls, feat_matrix=X, feat_df=FeatureDF):
    ref_idxs = []
    for idx, lbl in enumerate(ref_lbls):
        ref_idxs.append(feat_df.index[feat_df["Batch"] == lbl].tolist()[0])
    feature_values = list(feat_matrix[ref_idxs, :].mean(axis=0))
    return feature_values


print("--- cell ok ---")
def create_plot_object(
    eigenvector,
    feature_values,
    bar_color,
    feature_name,
    show_legend=True,
    feature_df=FeatureDF,
    N_PCs=3,
):
    # compute feature contribution
    contribution = [p * x for p, x in zip(eigenvector, feature_values)]

    # create plot object for the contribution
    plot_object = go.Bar(
        x=feature_df.columns.tolist()[2:],
        y=contribution,
        name=feature_name,
        marker=dict(color=bar_color),
        showlegend=show_legend,
    )
    return plot_object


print("--- create_plot_object defined ---")
def plot_feature_contributions(
    outlier_feature_vals, ref_feature_vals, eigenvector_matrix=P
):

    # initialize contribution plots
    fig = make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.03)

    # create plot objects for contributions in first 3 PCs
    for i_component in range(3):

        # flag, whether to show legend
        show_flag = True if i_component == 0 else False

        # eigenvector for PC
        eigenvector = eigenvector_matrix[i_component, :]

        # create plot object for outlier batch(es)
        outlier_obj = create_plot_object(
            eigenvector,
            outlier_feature_vals,
            "royalblue",
            "Outlier",
            show_legend=show_flag,
        )

        # create plot object for inspec batch(es)
        ref_obj = create_plot_object(
            eigenvector,
            ref_feature_vals,
            "lightsteelblue",
            "In-spec",
            show_legend=show_flag,
        )

        # create bar plot for PC contribution
        fig.append_trace(outlier_obj, i_component + 1, 1)
        fig.append_trace(ref_obj, i_component + 1, 1)

        # label y axis appropriately
        fig["layout"][f"yaxis{i_component+1}"].update(
            title=f"PC{i_component+1} Contribution"
        )

    # layout of bar plot
    fig["layout"].update(width=800, height=900, xaxis=dict(tickangle=45))

    # show bar plot
    fig.show()
    return

print("--- plot_feature_contributions defined ---")
############### <<<<< DO NOT CHANGE SECTION - END <<<<< #############

# call function extract feature values
outlier_feature_values = extract_feature_values(outlier_batches)
inspec_feature_values = extract_feature_values(inspec_batches)

# call function to plot contributions
plot_feature_contributions(outlier_feature_values, inspec_feature_values)<jupyter_output><empty_output><jupyter_text>**Highlight outlier rows in FeatureDF.csv**<jupyter_code># read the save feature dataframe
feature_df = pd.read_csv("FeatureDF.csv")
feature_df = feature_df.astype({"Batch": int})

# Filter on index+1 in BRANDS_OF_INTEREST (because brand in FeatureDF.csv is an integer)
brand_df = feature_df[feature_df["Brand"] == BRANDS_OF_INTEREST.index(brand) + 1]

# row turns red for batch flagged as outlier
def highlight_outliers(s):
    if s.Batch in outlier_batches:
        return ["background-color: red"] * 10
    else:
        return [""] * 10


brand_df.style.apply(highlight_outliers, axis=1)
from IPython.display import display, Markdown

display(Markdown("### QUESTIONS:"))
display(Markdown("**1. Why do you think there are bad batches?**"))
# =========== STUDENT BEGIN ==========
display(Markdown(
"""
@@@ Your answer here @@@
"""))
# =========== STUDENT END ==========

display(Markdown("**2. Does it seem like there is equipment failure at some point?**"))
# =========== STUDENT BEGIN ==========
display(Markdown(
"""
@@@ Your answer here @@@
"""))
# =========== STUDENT END ==========

display(Markdown("**3. What would you do to determine the root cause of a process or equipment failure?**"))
# =========== STUDENT BEGIN ==========
display(Markdown(
"""
@@@ Your answer here @@@
"""))
# =========== STUDENT END ==========<jupyter_output><empty_output><jupyter_text>--- Part 6. Apply clustering to Principal Components (BONUS)---There are many clustering algorithms which are available. One of the most commonly used is **k-means clustering**.  6a. Apply k-means clustering on the principal componentsWe use the cumulative variance plot to decide on an appropriate number ("most dominant") of principal components, and we perform k-means clustering only those components.<jupyter_code># base this on CPV plot
n_principal_components = 4

# number of clusters = number of available beer brands
K = len(FeatureDF["Brand"].unique())

# apply kmeans on principal components
kmeans = KMeans(n_clusters=K).fit(T[:, :n_principal_components])

# find which cluster each beer batch belongs to
cluster_id = kmeans.labels_<jupyter_output><empty_output><jupyter_text>6b. Plot principal components and color by kmeans cluster<jupyter_code># plot data for all clusters
obj_list = []

all_clusters = list(set(cluster_id))
for i_cluster, cluster in enumerate(all_clusters):

    boolean = cluster_id == cluster

    # create plotly object
    plot_obj = go.Scatter3d(
        x=T[boolean, 0],
        y=T[boolean, 1],
        z=T[boolean, 2],
        text=list(FeatureDF.loc[boolean, "Batch"]),
        hoverinfo="text",
        hovertemplate = 'Batch %{text:2.0}',
        name=f"Cluster {i_cluster+1}",
        mode="markers",
        marker=dict(size=7, opacity=0.75),
    )

    # append object to list
    obj_list.append(plot_obj)

layout = go.Layout(
    title="<br>K-Cluster Scatterplot",
    title_x=0.5,
    scene=dict(
        xaxis=dict(title="PC 1"), yaxis=dict(title="PC 2"), zaxis=dict(title="PC 3")
    ),
    showlegend=True,
    margin=dict(l=0, r=0, b=0, t=0),
    legend=dict(x=1.0, y=0.5),
)

fig = go.Figure(data=obj_list, layout=layout)
fig.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10608452081897606, 'fraction_numerical': 0.010226085410043795, 'mean_word_length': 3.904491931966856, 'pattern_counts': {'":': 2, '<': 114, '<?xml version=': 0, '>': 121, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '253892', 'license': 'Apache-2.0', 'n_tokens_mistral': 14019, 'n_tokens_neox': 13034, 'n_words': 4622, 'path': 'NB3_EXERCISE_Beer_PCA.ipynb', 'repo_name': 'academic-hub/brewery-exercises', 'size': '44983'}