Multiblock batch PLS: the FMC batch dryer#

Case study for issue #154.

An agricultural chemical is dried in an industrial batch dryer. Wet cake, the solid with its embedded solvent, is charged to the dryer and dried through three recipe phases: the solvent is collected in a side tank, the temperature is ramped, and the batch is cooled down. Chemical changes take place in the solid during drying, and the operators can adjust a few set points. Ten trajectories are recorded over each batch, the clock time at each aligned sample is carried as an eleventh, and three more blocks describe each batch with one row: the chemistry of the cake before the batch (Zchem, eleven measurements), the weight of the cake with eight landmarks of the batch’s own trajectories (Zop, nine values: the collector level and the dryer temperature at the end of the first phase, the peak temperature, the length of each phase and of the high-speed agitation, and the slope of the temperature ramp), and eight final quality attributes (Y). This is the multiblock case study of Garcia-Munoz and co-workers (2003).

The questions are the ones a plant asks in this order: what does product quality look like, do the initial conditions explain it, what do the trajectories add, and which batches deserve a closer look. The original course material answers them with a ladder of models, two components each, and this script follows the same ladder, then reads the block scores of the final model for the batches whose trajectories say off-specification while the product was on-specification.

The complete script is fmc_multiblock_batch_pls.py in this directory:

uv run python docs/user_guide/case_studies/batch/fmc_multiblock_batch_pls.py --output-dir case-study-output/fmc

It prints the numbers quoted below and writes every figure as an HTML file to the output directory.

Data#

Industrial batch dryer: a workbook with the four blocks over 59 batches. The trajectories were aligned within each of the three phases before the data were archived, to 325 samples per batch, and ClockTime, the wall time at each aligned sample, is carried along as an eleventh trajectory: after alignment it is no longer a clock but a record of how much each batch was stretched or compressed, which is itself information about the batch. (For raw, unaligned data see process_improve.batch.batch_dtw(); the unaligned trajectories of this same process are bundled as process_improve.batch.load_dryer().)

Thirteen batches have no chemistry measurements. The original study excluded them, and process_improve.batch.load_fmc() returns their identifiers as missing_chemistry so the exclusion can be reproduced. The remaining 46 batches still contain genuine missing cells: 19 in the quality block, one in the chemistry block, and 1220 in the trajectories of ten batches. The process_improve.multivariate estimators handle missing values through their NIPALS path, which is why this case study uses PCA, PLS and MBPLS directly; the batch classes BatchPCA and BatchPLS require complete data.

N_COMPONENTS = 2
CONF_LEVEL = 0.95
QUALITY_GROUP = [61, 14]  # one batch from each group in the quality score plot
OPERATING_OUTLIER = 20  # stands out on the operating conditions and on the trajectories
TRAJECTORY_BATCHES = [13, 5, 7]  # batches examined in the batch PLS
DISPOSITION = {"good": 33, "abnormal": 61, "high solvent": 71}  # the plant's classes: the last batch number of each
HIGHLIGHT = '{"color": "red", "width": 4}'  # Plotly line style, JSON-encoded, for the highlighted batches
NEIGHBOUR = '{"color": "teal", "width": 3}'  # the same, for the batches a highlighted batch is compared with
LABELS = {"show_labels": True}

def load_and_exclude(url: str | None = None) -> Bunch:
    """Download the four blocks and drop the batches without chemistry data."""
    fmc = load_fmc(url=url)
    keep = [batch_id for batch_id in fmc.batch_ids if batch_id not in fmc.missing_chemistry]
    data = Bunch(
        X={batch_id: fmc.X[batch_id] for batch_id in keep},
        Y=fmc.Y.loc[keep],
        Zop=fmc.Zop.loc[keep],
        Zchem=fmc.Zchem.loc[keep],
    )
    incomplete = [batch_id for batch_id, batch in data.X.items() if batch.isna().any().any()]
    print(
        f"{len(keep)} batches kept; missing cells: Y {int(data.Y.isna().sum().sum())}, Zchem {int(data.Zchem.isna().sum().sum())}, X in batches {incomplete}"
    )
    return data



46 batches kept; missing cells: Y 19, Zchem 1, X in batches [20, 22, 27, 28, 31, 55, 60, 61, 67, 71]

def plot_raw(trajectories: dict, tag: str, highlight: list[int]) -> go.Figure:
    """Overlay one tag for every batch, with the batches of interest drawn in red."""
    return plot_all_batches_per_tag(
        trajectories, tag, batches_to_highlight={HIGHLIGHT: highlight}, extra_info=f"highlighting {highlight}"
    )

Characterising product quality#

A PCA on the eight quality attributes shows how the batches group in quality space before any process data are involved.

def pca_on_quality(quality: pd.DataFrame) -> tuple[PCA, pd.DataFrame]:
    """Two-component PCA of the quality block, to see how the batches group in quality space."""
    y_scaled = MCUVScaler().fit_transform(quality)  # missing cells pass through; PCA switches to NIPALS
    model = PCA(n_components=N_COMPONENTS).fit(y_scaled)
    print(f"PCA on Y: R2 cumulative = {cumulative(model.r2_cumulative_)}")
    contributions = model.score_contributions(y_scaled, component=1)
    for batch_id in QUALITY_GROUP:
        print(f"batch {batch_id}: t1 contributions = " + describe(contributions.loc[batch_id]))
    return model, y_scaled


def describe(values: pd.Series) -> str:
    """Format a contribution vector, naming the cells that are missing in the data."""
    return ", ".join(f"{name} {'missing' if pd.isna(value) else f'{value:+.2f}'}" for name, value in values.items())


def cumulative(values: pd.Series | np.ndarray) -> str:
    """Format a cumulative R2 vector."""
    return ", ".join(f"{value:.3f}" for value in np.asarray(values, dtype=float))


def plot_bars(values: pd.Series, title: str) -> go.Figure:
    """Bar chart of one contribution vector."""
    fig = go.Figure(go.Bar(x=[str(name) for name in values.index], y=values.to_numpy(dtype=float)))
    fig.update_layout(title=title, yaxis_title="Contribution")
    return fig



PCA on Y: R2 cumulative = 0.500, 0.703
batch 61: t1 contributions = Y1 -0.76, Y2 missing, Y4 -1.04, Y6 -1.05, Y9 -0.10, Y10 -0.77, Y11 -1.07, SolventConc +0.02
batch 14: t1 contributions = Y1 +0.77, Y2 +0.28, Y4 +0.73, Y6 +0.78, Y9 -0.04, Y10 +0.84, Y11 +0.15, SolventConc -0.11

Two components explain 70% of the quality block and the score plot shows two groups of batches. Batches 61 and 14 are one member of each group; their \(t_1\) contributions are mirror images, with the same attributes (Y1, Y4, Y6, Y10 and Y11) low in one group and high in the other, so the first component is a general quality level rather than a trade-off between attributes. A missing quality cell simply has no contribution.

Effect of the initial conditions#

def pls_from_initial_conditions(data: Bunch, y_scaled: pd.DataFrame) -> tuple[PLS, PLS, pd.DataFrame]:
    """PLS from each initial-condition block to quality, one block at a time.

    The blocks are scaled with ``MCUVScaler`` first and the models fitted with
    ``scale=False``, so every later contribution plot works in the same scaled
    space as the model.
    """
    zchem_scaled = MCUVScaler().fit_transform(data.Zchem)
    zop_scaled = MCUVScaler().fit_transform(data.Zop)
    pls_chem = PLS(n_components=N_COMPONENTS, scale=False).fit(zchem_scaled, y_scaled)
    pls_op = PLS(n_components=N_COMPONENTS, scale=False).fit(zop_scaled, y_scaled)
    print(f"PLS Zchem -> Y: R2Y cumulative = {cumulative(pls_chem.r2_cumulative_)}")
    print(f"PLS Zop -> Y: R2Y cumulative = {cumulative(pls_op.r2_cumulative_)}")
    contributions = pls_op.score_contributions(zop_scaled, component=1).loc[OPERATING_OUTLIER]
    print(
        f"batch {OPERATING_OUTLIER} on Zop: t1 contributions = "
        + ", ".join(f"{name} {value:+.2f}" for name, value in contributions.items())
    )
    return pls_chem, pls_op, zop_scaled



PLS Zchem -> Y: R2Y cumulative = 0.163, 0.222
PLS Zop -> Y: R2Y cumulative = 0.207, 0.262
batch 20 on Zop: t1 contributions = Level1 -0.33, Temp1 +0.02, Temp2 -0.02, Time4 -1.01, Time1 -0.14,
                                    Time2 -1.42, Time3 -0.15, TempSlope -1.27, WgtCake +0.19

Each initial-condition block on its own explains about a quarter of the quality block. Batch 20 stands out in the operating-condition model through its recipe timings (Time2, Time4) and temperature slope, which is worth remembering when it turns up again in the trajectory models.

Multiblock PLS on the initial conditions#

The two blocks can be modelled together. MBPLS scales each block on its own and then weights it by \(1/\sqrt{K_b}\), so the eleven chemistry columns and the nine operating columns pull on the super-score with equal total weight; the block scores show what each block contributes, and the super-score plot shows the batches in the combined space.

def mbpls_on_initial_conditions(data: Bunch) -> MBPLS:
    """Multiblock PLS from both initial-condition blocks to quality.

    Each block is scaled on its own and then weighted by 1 / sqrt(K_b), so the
    eleven chemistry columns and the nine operating columns pull on the
    super-score with equal total weight.
    """
    blocks = {"Zchem": data.Zchem, "Zop": data.Zop}
    model = MBPLS(n_components=N_COMPONENTS).fit(blocks, data.Y)
    print(
        f"MBPLS Z -> Y: R2Y cumulative = {cumulative(model.r2_y_cumulative_)}; R2X per block after {N_COMPONENTS} components = "
        + ", ".join(f"{name} {value:.3f}" for name, value in model.r2_x_per_block_cumulative_.iloc[:, -1].items())
    )
    for name, block_contributions in model.score_contributions(blocks, component=1).items():
        row = block_contributions.loc[OPERATING_OUTLIER]
        print(
            f"batch {OPERATING_OUTLIER}, block {name}: t1 contributions = "
            + ", ".join(f"{col} {value:+.2f}" for col, value in row.items())
        )
    return model



MBPLS Z -> Y: R2Y cumulative = 0.292, 0.364; R2X per block after 2 components = Zchem 0.296, Zop 0.356
batch 20, block Zchem: t1 contributions = Z1 -0.08, Z2 +0.03, Z3 -0.00, Z4 +0.04, Z5 +0.00, Z6 +0.03, Z7 -0.05, ...
batch 20, block Zop: t1 contributions = Level1 -0.08, Temp1 +0.02, Temp2 -0.01, Time4 -0.25, Time1 -0.04, Time2 -0.38, ...

Together the blocks explain 36% of the quality block, more than either alone, and the per-block contributions show that batch 20 is unusual in its operating conditions, not in its chemistry.

The trajectories alone#

The trajectories are unfolded batchwise with process_improve.batch.dict_to_wide(): one row per batch of 11 tags (the ten process measurements and ClockTime) times 325 samples, 3575 columns. MCUVScaler returns flat column labels, so the 2-level (tag, sequence) index is re-attached after scaling; the batch plots read it.

def unfold_trajectories(trajectories: dict) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Unfold the trajectories batchwise and scale every column.

    One row per batch, 11 tags x 325 samples = 3575 columns: the ten process
    measurements and ``ClockTime``, the eleventh trajectory that records how
    each batch used time. ``MCUVScaler`` returns flat column labels, so the
    2-level ``(tag, sequence)`` index is re-attached; the batch plots need it.
    """
    wide = dict_to_wide(trajectories)
    x_scaled = MCUVScaler().fit_transform(wide)
    x_scaled.columns = wide.columns
    print(
        f"unfolded trajectories: {wide.shape[0]} batches x {wide.shape[1]} columns, {int(wide.isna().sum().sum())} missing cells"
    )
    return wide, x_scaled



def batch_pca_on_trajectories(x_scaled: pd.DataFrame) -> tuple[PCA, pd.DataFrame]:
    """Two-component batch PCA on the trajectories alone.

    Returns the model and the share of every cell in each batch's SPE. A batch
    with missing cells has its scores estimated from the observed cells, so its
    contributions are defined everywhere except at the cells it lacks.
    """
    model = PCA(n_components=N_COMPONENTS).fit(x_scaled)  # NIPALS, because of the missing cells
    print(f"batch PCA on X: R2 cumulative = {cumulative(model.r2_cumulative_)}")
    squared = model.spe_contributions(x_scaled) ** 2  # NaN only at the missing cells
    spe_share = squared.div(squared.sum(axis=1), axis=0) * 100
    by_tag = spe_share.loc[OPERATING_OUTLIER].groupby(level="tag", sort=False).sum()
    print(
        f"batch {OPERATING_OUTLIER}, above both limits: share of the SPE per tag = "
        + ", ".join(f"{tag} {share:.0f}%" for tag, share in by_tag.items())
    )
    return model, spe_share



unfolded trajectories: 46 batches x 3575 columns, 1340 missing cells
batch PCA on X: R2 cumulative = 0.231, 0.376
batch 20, above both limits: share of the SPE per tag = Agitator 3%, CTankLvl 4%, ClockTime 3%, D-Temp 9%,
    D-Temp-SP 6%, DiffPres 4%, DryPress 49%, J-Temp 6%, J-Temp-SP 3%, Power 6%, Torque 6%

Two components describe 38% of the batch-to-batch variation in the trajectories, and the time-varying loading plot shows where in the batch each component acts. Batch 20 is above both limits, and its record has gaps between samples 34 and 109. Its scores are estimated from the observed cells, so its SPE contributions are defined at every observed cell and absent at the missing ones, which carry no residual: half of its SPE sits in the dryer pressure, most of that in the first phase. The raw overlays of the dryer temperature, power and torque show the same batch against the rest.

Trajectories to quality#

def batch_pls_to_quality(x_scaled: pd.DataFrame, y_scaled: pd.DataFrame) -> tuple[PLS, pd.DataFrame]:
    """Two-component batch PLS from the unfolded trajectories to quality."""
    model = PLS(n_components=N_COMPONENTS, scale=False).fit(x_scaled, y_scaled)
    print(f"batch PLS X -> Y: R2Y cumulative = {cumulative(model.r2_cumulative_)}")
    contributions = model.score_contributions(x_scaled, component=1)
    by_tag = contributions.loc[TRAJECTORY_BATCHES[0]].groupby(level="tag", sort=False).sum()
    print(
        f"batch {TRAJECTORY_BATCHES[0]}: t1 contributions per tag = "
        + ", ".join(f"{tag} {value:+.1f}" for tag, value in by_tag.items())
    )
    return model, contributions



batch PLS X -> Y: R2Y cumulative = 0.266, 0.410
batch 13: t1 contributions per tag = Agitator -1.7, CTankLvl -8.0, ClockTime -8.1, D-Temp -4.7, D-Temp-SP -1.7,
          DiffPres -1.4, DryPress -1.0, J-Temp -1.1, J-Temp-SP -4.2, Power -2.9, Torque -2.7

The trajectories explain 41% of the quality block, more than the initial conditions did. Batch 13 is at one end of \(t_1\), and its contributions are spread over the tags with the clock time and the collector tank level (the pace of the batch and the solvent removal) leading; batches 5 and 7 are examined the same way. The observed-versus-predicted plot of SolventConc shows how well the residual solvent, the attribute the plant cares most about, follows from the trajectories.

Batch multiblock PLS#

The final model joins all three X blocks. The trajectory block enters as 3575 columns, and its \(1/\sqrt{K_b}\) weight keeps it from drowning out the two small blocks.

def batch_mbpls(data: Bunch, wide: pd.DataFrame) -> tuple[MBPLS, dict]:
    """Batch multiblock PLS: chemistry, operating conditions and trajectories to quality."""
    blocks = {"Zchem": data.Zchem, "Zop": data.Zop, "X": wide}
    model = MBPLS(n_components=N_COMPONENTS).fit(blocks, data.Y)
    print(
        f"batch MBPLS: R2Y cumulative = {cumulative(model.r2_y_cumulative_)}; R2X per block after {N_COMPONENTS} components = "
        + ", ".join(f"{name} {value:.3f}" for name, value in model.r2_x_per_block_cumulative_.iloc[:, -1].items())
    )
    print("super VIP per block: " + ", ".join(f"{name} {value:.2f}" for name, value in model.super_vip_.items()))
    return model, blocks



batch MBPLS: R2Y cumulative = 0.370, 0.472; R2X per block after 2 components = Zchem 0.233, Zop 0.304, X 0.259
super VIP per block: Zchem 0.86, Zop 1.07, X 1.06

The combined model explains 47% of the quality block, and the super VIP puts the operating conditions and the trajectories about level and the chemistry last. This is the model to build the stagewise monitoring and the final-quality prediction on: the super-score plot places every batch in one space, the block scores say whether a batch is unusual in its chemistry, its operation or its trajectories, and the X-block contributions of a batch, drawn with process_improve.batch.unfolded_contribution_plot(), name the tags and the phase.

Off-specification trajectories, on-specification product#

The block scores are read next. Each batch is placed, block by block, with the group (good or abnormal, by the plant’s disposition) whose average point is nearer in that block’s score plot. Four batches classed good are placed with the abnormal batches by the trajectory block and with the good batches by both initial-condition blocks; batch 5, placed with the abnormal batches by the operating-condition block as well, is left aside. The four are compared with their nearest abnormal neighbours in the trajectory block through the contribution from the neighbours’ average to theirs in the operating-condition block.

def disposition(batch_ids: list) -> pd.Series:
    """Return the plant's disposition of each batch, which is encoded in the batch numbering."""
    edges = [0, *DISPOSITION.values()]
    return pd.Series(pd.cut(batch_ids, bins=edges, labels=list(DISPOSITION)).astype(str), index=batch_ids)


def nearer_group(scores: pd.DataFrame, groups: pd.Series) -> pd.Series:
    """Place each batch with the group, good or abnormal, whose average point is nearer in this score plot."""
    centres = {name: scores.loc[groups == name].mean() for name in ("good", "abnormal")}
    return pd.DataFrame({name: ((scores - centre) ** 2).sum(axis=1) for name, centre in centres.items()}).idxmin(axis=1)


def off_spec_trajectories_on_spec_product(model: MBPLS, blocks: dict) -> Bunch:
    """Find the batches classed good whose trajectories sit with the abnormal batches, and what set them apart.

    Every block of the batch multiblock PLS has its own score plot, and a
    batch is placed in each of them with the group whose average point is
    nearer. The batches classed good that the trajectory block places with
    the abnormal batches, while both initial-condition blocks place them with
    the good ones, are compared with their nearest abnormal neighbours in the
    trajectory block: the contribution from the neighbours' average to theirs
    in the operating-condition block names what differed.
    """
    groups = disposition(list(model.super_scores_.index))
    placed = pd.DataFrame({name: nearer_group(scores, groups) for name, scores in model.block_scores_.items()})
    with_abnormal = [b for b in placed.index if groups[b] == "good" and placed.loc[b, "X"] == "abnormal"]
    anomalous = [b for b in with_abnormal if (placed.loc[b, ["Zchem", "Zop"]] == "good").all()]
    print(f"batches classed good that the trajectory block places with the abnormal batches: {with_abnormal}")
    print(f"of these, placed with the good batches by both initial-condition blocks: {anomalous}")
    x_scores = model.block_scores_["X"]
    abnormal = x_scores.loc[groups == "abnormal"]
    nearest = {int(b) for a in anomalous for b in ((abnormal - x_scores.loc[a]) ** 2).sum(axis=1).nsmallest(2).index}
    neighbours = sorted(nearest)
    contributions = model.score_contributions(blocks, component=1)["Zop"]
    move = contributions.loc[anomalous].mean() - contributions.loc[neighbours].mean()
    print(f"their nearest abnormal batches in the trajectory block: {neighbours}")
    print(
        "Zop contribution from the neighbours' average to the anomalous batches' average: "
        + ", ".join(f"{name} {value:+.2f}" for name, value in move.items())
    )
    return Bunch(anomalous=anomalous, neighbours=neighbours, placed=placed, zop_move=move)


def plot_block_scores(model: MBPLS, mark: list[int]) -> go.Figure:
    """Draw the score plot of every block side by side, coloured by the plant's disposition, with `mark` labelled."""
    groups = disposition(list(model.super_scores_.index))
    colours = {"good": "steelblue", "abnormal": "purple", "high solvent": "teal"}
    fig = make_subplots(
        rows=1, cols=len(model.block_scores_), subplot_titles=[f"{n} block" for n in model.block_scores_]
    )
    for col, (name, scores) in enumerate(model.block_scores_.items(), start=1):
        for label, colour in colours.items():
            members = [b for b in scores.index if groups[b] == label and b not in mark]
            fig.add_trace(
                go.Scatter(
                    x=scores.loc[members].iloc[:, 0],
                    y=scores.loc[members].iloc[:, 1],
                    mode="markers",
                    name=f"classed {label}",
                    marker={"color": colour},
                    text=members,
                    hovertemplate="batch %{text}",
                    showlegend=col == 1,
                ),
                row=1,
                col=col,
            )
        fig.add_trace(
            go.Scatter(
                x=scores.loc[mark].iloc[:, 0],
                y=scores.loc[mark].iloc[:, 1],
                mode="markers+text",
                text=[str(b) for b in mark],
                textposition="top right",
                name="anomalous",
                marker={"color": "red", "size": 10},
                showlegend=col == 1,
            ),
            row=1,
            col=col,
        )
        r2 = np.diff([0.0, *model.r2_x_per_block_cumulative_.loc[name].to_numpy(dtype=float)])
        fig.update_xaxes(title_text=f"t1 [{r2[0]:.1%}]", row=1, col=col)
        fig.update_yaxes(title_text=f"t2 [{r2[1]:.1%}]", row=1, col=col)
    fig.update_layout(title="Block scores of the batch multiblock PLS", height=420)
    return fig



batches classed good that the trajectory block places with the abnormal batches: [2, 3, 5, 6, 7]
of these, placed with the good batches by both initial-condition blocks: [2, 3, 6, 7]
their nearest abnormal batches in the trajectory block: [42, 43, 44, 47, 50]
Zop contribution from the neighbours' average to the anomalous batches' average: Level1 -0.05, Temp1 +0.02,
    Temp2 +0.00, Time4 +0.05, Time1 -0.00, Time2 +0.06, Time3 +0.11, TempSlope +0.06, WgtCake -0.05

The four batches share their neighbours’ heavy charge and high collector level (WgtCake and Level1 pull towards the abnormal side) and differ from them in the later phases: a longer cool-down (Time3, the largest single contribution), a shorter and steeper temperature ramp (Time2, TempSlope) and a shorter high-speed agitation (Time4). The peak temperature set point does not differ, so this is not a set point that was moved but a difference in how long each phase was run. Whether the later phases were run that way to correct for the first, the record does not say, and the model describes how the batches co-varied, not cause and effect (Nomikos and MacGregor, 1995). Garcia-Munoz (2004) reads the same four batches the same way in Appendix 1 of the thesis.

Where to go next#

The course notes suggest replacing the raw trajectories with feature blocks (timings, temperatures, impeller and pressure summaries) so that the model can be read by phase, and building the online monitoring model with process_improve.batch.BatchMonitor once the missing cells have been filled in or the incomplete batches removed.

Running the script#

    """Run the whole case study and write its figures."""
    parser = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
    parser.add_argument("--output-dir", type=pathlib.Path, default=pathlib.Path("case-study-output/fmc"))
    parser.add_argument("--data-url", default=None, help="override the openmv.net URL (for example a file:// copy)")
    args = parser.parse_args(argv)
    args.output_dir.mkdir(parents=True, exist_ok=True)

    def save(fig: go.Figure, name: str) -> None:
        fig.write_html(args.output_dir / f"{name}.html", include_plotlyjs="cdn")

    data = load_and_exclude(args.data_url)
    for tag in ("D-Temp", "J-Temp", "CTankLvl", "ClockTime"):
        save(plot_raw(data.X, tag, [OPERATING_OUTLIER]), f"raw-{tag}")

    pca_y, y_scaled = pca_on_quality(data.Y)
    save(pca_y.score_plot(settings=LABELS), "quality-pca-scores")
    y_contributions = pca_y.score_contributions(y_scaled, component=1)
    for batch_id in QUALITY_GROUP:
        save(
            plot_bars(y_contributions.loc[batch_id], f"t1 contributions of batch {batch_id} (PCA on Y)"),
            f"quality-pca-contributions-{batch_id}",
        )

    _pls_chem, pls_op, zop_scaled = pls_from_initial_conditions(data, y_scaled)
    save(pls_op.score_plot(settings=LABELS), "pls-zop-scores")
    save(
        plot_bars(
            pls_op.score_contributions(zop_scaled, component=1).loc[OPERATING_OUTLIER],
            f"t1 contributions of batch {OPERATING_OUTLIER} (PLS Zop)",
        ),
        "pls-zop-contributions-20",
    )

    mbpls_z = mbpls_on_initial_conditions(data)
    save(mbpls_z.super_score_plot(), "mbpls-z-super-scores")
    save(mbpls_z.super_weights_bar_plot(component=1), "mbpls-z-super-weights")

    wide, x_scaled = unfold_trajectories(data.X)
    pca_x, spe_share = batch_pca_on_trajectories(x_scaled)
    save(pca_x.score_plot(settings=LABELS), "batch-pca-scores")
    save(pca_x.spe_plot(settings=LABELS), "batch-pca-spe")
    save(time_varying_loading_plot(pca_x, component=1), "batch-pca-loadings-p1")
    save(
        unfolded_contribution_plot(spe_share.fillna(0.0), OPERATING_OUTLIER, by_tag=True),
        f"batch-pca-spe-contributions-{OPERATING_OUTLIER}",
    )
    for tag in ("D-Temp", "Power", "Torque"):
        save(plot_raw(data.X, tag, [OPERATING_OUTLIER]), f"raw-{tag}-batch-20")

    pls_x, x_contributions = batch_pls_to_quality(x_scaled, y_scaled)
    save(pls_x.score_plot(settings=LABELS), "batch-pls-scores")
    save(unfolded_contribution_plot(x_contributions, TRAJECTORY_BATCHES[0]), "batch-pls-contributions-13")
    for batch_id in TRAJECTORY_BATCHES:
        save(plot_raw(data.X, "D-Temp", [batch_id]), f"raw-D-Temp-batch-{batch_id}")
    save(
        pls_x.predictions_vs_observed_plot(y_observed=y_scaled, variable="SolventConc"),
        "batch-pls-observed-vs-predicted",
    )

    mbpls_x, blocks = batch_mbpls(data, wide)
    save(mbpls_x.super_score_plot(), "batch-mbpls-super-scores")
    save(mbpls_x.super_weights_bar_plot(component=1), "batch-mbpls-super-weights")
    save(
        unfolded_contribution_plot(mbpls_x.score_contributions(blocks, component=1)["X"], TRAJECTORY_BATCHES[0]),
        "batch-mbpls-x-contributions-13",
    )
    save(mbpls_x.predictions_vs_observed_plot(data.Y, variable="SolventConc"), "batch-mbpls-observed-vs-predicted")

    found = off_spec_trajectories_on_spec_product(mbpls_x, blocks)
    save(plot_block_scores(mbpls_x, found.anomalous), "batch-mbpls-block-scores")
    save(
        plot_bars(found.zop_move, "Zop contribution from the neighbours' average to the anomalous batches'"),
        "batch-mbpls-zop-move",
    )
    for tag in ("CTankLvl", "ClockTime", "D-Temp", "D-Temp-SP"):
        fig = plot_all_batches_per_tag(
            data.X,
            tag,
            batches_to_highlight={HIGHLIGHT: found.anomalous, NEIGHBOUR: found.neighbours},
            extra_info=f"anomalous {found.anomalous} in red, their neighbours {found.neighbours} in teal",
        )
        save(fig, f"raw-{tag}-anomalous")
    print(f"figures written to {args.output_dir}")
    return 0

References#