This vignette introduces the main functions of the
PEAXAI package.
PEAXAI provides a methodology that integrates Data
Envelopment Analysis (DEA) with Machine Learning (ML) algorithms to
train a classifier that estimates the probability that each
decision-making unit (DMU) is (in)efficient.
Once the model is trained, the package offers several functions for post-hoc analysis, including:
These tools facilitate benchmarking and policy simulations in a transparent and explainable framework.
To illustrate the functionality of the PEAXAI package,
we use a dataset of 917 companies in the food industry sector operating
in Spain. The data were retrieved from the SABI (Sistema de Análisis de
Balances Ibéricos) database for the year 2023 and include only firms
with more than 50 employees.
Each row corresponds to a single firm and includes both financial and operational variables relevant for productivity and efficiency analysis.
The output variable used to measure performance is:
operating_income: Operating income in millions of
euros.The input variables include:
total_assets: Total assets (in millions of euros).employees: Number of employees.fixed_assets: Tangible fixed assets (in millions of
euros).personnel_expenses: Personnel-related costs (in
millions of euros).Additionally, the variable autonomous_community
indicates the geographical location of each firm within one of Spain’s
17 autonomous communities and 2 autonomous cities. This allows the
analysis to reflect regional institutional and market heterogeneity.
The dataset exhibits a wide dispersion across firms, including both medium-sized and large enterprises. This heterogeneity poses a realistic challenge for evaluating efficiency and provides a rich testing ground for model explainability.
For illustration purposes, this vignette focuses on a subset of the
data: the 97 firms located in the Autonomous Community of Valencia
(Comunidad Valenciana). We drop the autonomous_community
column, so that the remaining columns are exactly the four inputs and
the output.
data <- subset(
firms,
autonomous_community == "Comunidad Valenciana",
select = -autonomous_community
)
rm(firms)Quick structure and summary.
str(data)
#> 'data.frame': 97 obs. of 5 variables:
#> $ total_assets : num 259 185 129 104 153 ...
#> $ employees : num 70 261 354 968 1076 ...
#> $ fixed_assets : num 84.15 21.75 5.76 30.55 103.76 ...
#> $ personnel_expenses: num 16.7 19.9 13.7 36.8 31.8 ...
#> $ operating_income : num 461 358 357 294 268 ...
summary(data)
#> total_assets employees fixed_assets personnel_expenses
#> Min. : 1.537 Min. : 50.0 Min. : 0.1421 Min. : 1.037
#> 1st Qu.: 8.989 1st Qu.: 75.0 1st Qu.: 2.6797 1st Qu.: 2.167
#> Median : 24.555 Median : 98.0 Median : 6.2578 Median : 3.059
#> Mean : 41.030 Mean : 201.1 Mean : 15.2799 Mean : 6.757
#> 3rd Qu.: 52.409 3rd Qu.: 240.0 3rd Qu.: 20.0960 3rd Qu.: 8.244
#> Max. :258.825 Max. :1076.0 Max. :140.6890 Max. :36.789
#> operating_income
#> Min. : 2.382
#> 1st Qu.: 12.994
#> Median : 29.138
#> Mean : 62.307
#> 3rd Qu.: 72.688
#> Max. :460.578Determine inputs (x) and outputs (y) indices.
Determine the technology assumptions.
To address class imbalance, we apply SMOTE based on the best-practice frontier (assuming a convex production technology). First, we identify the combinations of DMUs that define the frontier. Then, we generate one of two types of synthetic DMUs from these combinations: efficient units (on the frontier) and near-frontier inefficient units (inside the frontier’s envelope).
When the efficient class is the minority, efficient-class SMOTE on the best-practice frontier helps the classifier learn the frontier more accurately. Optionally, if the observed imbalance still falls short of the target ratio, we also generate near-frontier inefficient synthetic units derived from the frontier geometry to sharpen the separation between efficient and inefficient regions.
Define machine learning methods and their tuning parameters. Here we configure a neural network (nnet) with a grid of hyperparameters
methods <- list(
"nnet" = list(
tuneGrid = expand.grid(
size = c(1, 5, 10, 20),
decay = 10^seq(-5, -1, by = 1)
),
preProcess = c("center", "scale"),
# --- arguments passed on to nnet ---
skip = TRUE,
maxit = 100,
MaxNWts = 100000,
trace = FALSE
)
)Parameters for controlling model training (k-fold cross-validation)
Classification metrics to optimize during training (priority order in case of ties)
Define hold-out sample to evaluate model performance on unseen data.
These samples are excluded from the cross-validation and hyperparameter
tuning process. By default, is NULL.
Set random seed for reproducibility (ensures consistent results across runs).
We now apply the full PEAXAI pipeline:
models <- PEAXAI_fitting(
data = data,
x = x,
y = y,
RTS = RTS,
imbalance_rate = imbalance_rate,
methods = methods,
trControl = trControl,
metric_priority = metric_priority,
hold_out = hold_out,
verbose = TRUE,
seed = seed
)
#> [1] "Computing maximal friends with 1 DMUs (step 1 of 15 )"
#> [1] "Computing maximal friends with 2 DMUs (step 2 of 15 )"
#> [1] "Computing maximal friends with 3 DMUs (step 3 of 15 )"
#> [1] "Computing maximal friends with 4 DMUs (step 4 of 15 )"
#> [1] "Computing maximal friends with 5 DMUs (step 5 of 15 )"
#> [1] "Computing maximal friends with 6 DMUs (step 6 of 15 )"
#> [1] "Computing maximal friends with 7 DMUs (step 7 of 15 )"
#> [1] "Computing maximal friends with 8 DMUs (step 8 of 15 )"
#> [1] "Computing maximal friends with 9 DMUs (step 9 of 15 )"
#> [1] "Computing maximal friends with 10 DMUs (step 10 of 15 )"
#> [1] "Computing maximal friends with 11 DMUs (step 11 of 15 )"
#> [1] "Computing maximal friends with 12 DMUs (step 12 of 15 )"
#> [1] "Computing maximal friends with 13 DMUs (step 13 of 15 )"
#> [1] "Computing maximal friends with 14 DMUs (step 14 of 15 )"
#> [1] "Computing maximal friends with 15 DMUs (step 15 of 15 )"
#> [1] "Computing maximal friends with 1 DMUs (step 1 of 14 )"
#> [1] "Computing maximal friends with 2 DMUs (step 2 of 14 )"
#> [1] "Computing maximal friends with 3 DMUs (step 3 of 14 )"
#> [1] "Computing maximal friends with 4 DMUs (step 4 of 14 )"
#> [1] "Computing maximal friends with 5 DMUs (step 5 of 14 )"
#> [1] "Computing maximal friends with 6 DMUs (step 6 of 14 )"
#> [1] "Computing maximal friends with 7 DMUs (step 7 of 14 )"
#> [1] "Computing maximal friends with 8 DMUs (step 8 of 14 )"
#> [1] "Computing maximal friends with 9 DMUs (step 9 of 14 )"
#> [1] "Computing maximal friends with 10 DMUs (step 10 of 14 )"
#> [1] "Computing maximal friends with 11 DMUs (step 11 of 14 )"
#> [1] "Computing maximal friends with 12 DMUs (step 12 of 14 )"
#> [1] "Computing maximal friends with 13 DMUs (step 13 of 14 )"
#> [1] "Computing maximal friends with 14 DMUs (step 14 of 14 )"
#> [1] "Computing maximal friends with 1 DMUs (step 1 of 12 )"
#> [1] "Computing maximal friends with 2 DMUs (step 2 of 12 )"
#> [1] "Computing maximal friends with 3 DMUs (step 3 of 12 )"
#> [1] "Computing maximal friends with 4 DMUs (step 4 of 12 )"
#> [1] "Computing maximal friends with 5 DMUs (step 5 of 12 )"
#> [1] "Computing maximal friends with 6 DMUs (step 6 of 12 )"
#> [1] "Computing maximal friends with 7 DMUs (step 7 of 12 )"
#> [1] "Computing maximal friends with 8 DMUs (step 8 of 12 )"
#> [1] "Computing maximal friends with 9 DMUs (step 9 of 12 )"
#> [1] "Computing maximal friends with 10 DMUs (step 10 of 12 )"
#> [1] "Computing maximal friends with 11 DMUs (step 11 of 12 )"
#> [1] "Computing maximal friends with 12 DMUs (step 12 of 12 )"
#> [1] "Computing maximal friends with 1 DMUs (step 1 of 14 )"
#> [1] "Computing maximal friends with 2 DMUs (step 2 of 14 )"
#> [1] "Computing maximal friends with 3 DMUs (step 3 of 14 )"
#> [1] "Computing maximal friends with 4 DMUs (step 4 of 14 )"
#> [1] "Computing maximal friends with 5 DMUs (step 5 of 14 )"
#> [1] "Computing maximal friends with 6 DMUs (step 6 of 14 )"
#> [1] "Computing maximal friends with 7 DMUs (step 7 of 14 )"
#> [1] "Computing maximal friends with 8 DMUs (step 8 of 14 )"
#> [1] "Computing maximal friends with 9 DMUs (step 9 of 14 )"
#> [1] "Computing maximal friends with 10 DMUs (step 10 of 14 )"
#> [1] "Computing maximal friends with 11 DMUs (step 11 of 14 )"
#> [1] "Computing maximal friends with 12 DMUs (step 12 of 14 )"
#> [1] "Computing maximal friends with 13 DMUs (step 13 of 14 )"
#> [1] "Computing maximal friends with 14 DMUs (step 14 of 14 )"
#> [1] "Computing maximal friends with 1 DMUs (step 1 of 13 )"
#> [1] "Computing maximal friends with 2 DMUs (step 2 of 13 )"
#> [1] "Computing maximal friends with 3 DMUs (step 3 of 13 )"
#> [1] "Computing maximal friends with 4 DMUs (step 4 of 13 )"
#> [1] "Computing maximal friends with 5 DMUs (step 5 of 13 )"
#> [1] "Computing maximal friends with 6 DMUs (step 6 of 13 )"
#> [1] "Computing maximal friends with 7 DMUs (step 7 of 13 )"
#> [1] "Computing maximal friends with 8 DMUs (step 8 of 13 )"
#> [1] "Computing maximal friends with 9 DMUs (step 9 of 13 )"
#> [1] "Computing maximal friends with 10 DMUs (step 10 of 13 )"
#> [1] "Computing maximal friends with 11 DMUs (step 11 of 13 )"
#> [1] "Computing maximal friends with 12 DMUs (step 12 of 13 )"
#> [1] "Computing maximal friends with 13 DMUs (step 13 of 13 )"
#> [1] "Computing maximal friends with 1 DMUs (step 1 of 14 )"
#> [1] "Computing maximal friends with 2 DMUs (step 2 of 14 )"
#> [1] "Computing maximal friends with 3 DMUs (step 3 of 14 )"
#> [1] "Computing maximal friends with 4 DMUs (step 4 of 14 )"
#> [1] "Computing maximal friends with 5 DMUs (step 5 of 14 )"
#> [1] "Computing maximal friends with 6 DMUs (step 6 of 14 )"
#> [1] "Computing maximal friends with 7 DMUs (step 7 of 14 )"
#> [1] "Computing maximal friends with 8 DMUs (step 8 of 14 )"
#> [1] "Computing maximal friends with 9 DMUs (step 9 of 14 )"
#> [1] "Computing maximal friends with 10 DMUs (step 10 of 14 )"
#> [1] "Computing maximal friends with 11 DMUs (step 11 of 14 )"
#> [1] "Computing maximal friends with 12 DMUs (step 12 of 14 )"
#> [1] "Computing maximal friends with 13 DMUs (step 13 of 14 )"
#> [1] "Computing maximal friends with 14 DMUs (step 14 of 14 )"#> $nnet
#> Neural Network
#>
#> 127 samples
#> 5 predictor
#> 2 classes: 'efficient', 'not_efficient'
#>
#> Pre-processing: centered (5), scaled (5)
#> Resampling: None#> $nnet
#> Imbalance_rate size decay Accuracy Kappa Recall Specificity Precision F1
#> 141 0.35 1 1e-05 0.94 0.74 0.8 0.96 0.85 0.77
#> Balanced_Accuracy G_mean ROC_AUC PR_AUC Cross_Entropy
#> 141 0.88 0.86 0.98 0.66 0.25
#> Cross_Entropy_Efficient_class Cross_Entropy_not_Efficient_class AccuracySD
#> 141 0.3 0.24 0.02
#> KappaSD RecallSD SpecificitySD PrecisionSD F1SD Balanced_AccuracySD
#> 141 0.16 0.3 0.03 0.14 0.16 0.13
#> G_meanSD ROC_AUCSD PR_AUCSD Cross_EntropySD Cross_Entropy_Efficient_classSD
#> 141 0.17 0.02 0.01 0.17 0.42
#> Cross_Entropy_not_Efficient_classSD diff_imbalance
#> 141 0.25 0.8046This example uses SHAP (SHapley Additive exPlanations) computed with
the kernelshap package to estimate global feature importance with
PEAXAI_global_importance. We approximate Shapley values via
a background sample (bg_n) and obtain a global importance
score by aggregating absolute SHAP values across observations (and then
normalizing to relative importance). When bg_n is not
specified, a background sample of 200 observations is used. For
alternative XAI methods and options, see the package documentation.
Two datasets are supplied.
explain_data holds the DMUs on which the explanations
are computed. Passing the observed sample, as we do here, yields
importances that represent the underlying problem as it is observed in
practice.
reference_data is the distribution the method perturbs
against, that is, the background sample for SHAP or the reference
distribution for a surrogate model. Passing
final_model$trainingData reflects what the model actually
learned during fitting, including the synthetic DMUs generated by SMOTE;
passing the observed sample restricts the background to real DMUs.
Both datasets must carry the variables the model was trained on under their exact names.
relative_importance <- PEAXAI_global_importance(
final_model = models[["best_model_fit"]][["nnet"]],
x = x,
y = y,
explain_data = data,
reference_data = data,
importance_method = importance_method
)
#> | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |========= | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |============ | 16% | |============ | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |=================================== | 49% | |=================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================== | 82% | |========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%#> total_assets employees fixed_assets personnel_expenses
#> importance 0.1297867 0.1702557 0.1298721 0.1053863
#> operating_income
#> importance 0.4646992We evaluate results at three efficiency cutoffs: 0.75, 0.85, and 0.95.
We define a global, model-aware direction for improving inputs/outputs.
relative_importance is the direction along which the
counterfactual analysis moves each DMU. It can be supplied by the user
or taken from PEAXAI_global_importance(), which returns a
single row and therefore applies one direction to every DMU. A matrix
with one row per DMU, such as the one returned by
PEAXAI_local_importance(), applies a DMU-specific direction
instead.
baseline sets the unit of that direction:
"mean", "median", "self" or
"ones".
For example, if the directional vector is custom: If it not possible
(or we do not want) to change a specific input (like
employees), we need to write:
relative_importance_custom <- t(matrix(
data = c(0.2, 0, 0.2, 0.2, 0.4),
))
relative_importance_custom <- as.data.frame(relative_importance_custom)
names(relative_importance_custom) <- names(data)[c(x,y)]#> total_assets employees fixed_assets personnel_expenses operating_income
#> 1 0.2 0 0.2 0.2 0.4But, if we do not know which direction define, we can use the relative importace by the model fitted:
Given the thresholds and directional vector, we compute feasible improvement targets for each DMU. Key controls:
n_expand = 0.25: expansion factor for the search
neighborhood.
n_grid = 50: grid resolution for exploring
adjustments.
max_y = 1, min_x = 1: bounds that cap
output expansion and input contraction.
targets <- PEAXAI_counterfactuals(
data = data,
x = x,
y = y,
final_model = models[["best_model_fit"]][["nnet"]],
efficiency_thresholds = efficiency_thresholds,
directional_vector = directional_vector,
n_expand = 0.25,
n_grid = 50,
max_y = 1,
min_x = 1
)head(targets[["0.85"]][["counterfactual_dataset"]], 10)
#> total_assets employees fixed_assets personnel_expenses operating_income
#> 1 258.53448 70.0000 84.15178 16.70692 460.5779
#> 2 184.97800 261.0000 21.75500 19.92400 358.1300
#> 3 128.93344 354.0000 5.75535 13.70937 356.9468
#> 4 102.95008 959.8593 30.07690 36.62022 300.8164
#> 5 141.33833 999.5421 99.32907 30.16345 332.9188
#> 6 155.56878 721.9301 76.24771 33.52334 287.0228
#> 7 108.60426 702.9214 50.13356 23.10993 254.5889
#> 8 237.83874 511.0362 132.86831 25.06555 303.4787
#> 9 53.41076 442.0000 24.98992 19.24766 170.1763
#> 10 81.40315 633.3956 19.34343 20.36560 224.0352We compute efficiency rankings at each threshold using two bases:
Predicted: ranks by the model’s predicted probability of
being efficient (rank_basis = predicted).
Attainable: ranks by attainable improvements/targets
implied by the model (rank_basis =
attainable).
ranking <- PEAXAI_ranking(
data = data,
x = x,
y = y,
final_model = models[["best_model_fit"]][["nnet"]],
rank_basis = "predicted"
)
ranking2 <- PEAXAI_ranking(
data = data,
x = x,
y = y,
final_model = models[["best_model_fit"]][["nnet"]],
efficiency_thresholds = efficiency_thresholds,
targets = targets,
rank_basis = "attainable"
)head(round(ranking, 4), 50)
#> Ranking DMU Probability_predicted
#> 1 1 1 1.0000
#> 2 2 2 1.0000
#> 3 3 3 1.0000
#> 4 4 17 1.0000
#> 5 5 20 1.0000
#> 6 6 46 1.0000
#> 7 7 9 0.9996
#> 8 8 62 0.9946
#> 9 9 92 0.9768
#> 10 10 18 0.8956
#> 11 11 36 0.8672
#> 12 12 56 0.7489
#> 13 13 75 0.7073
#> 14 14 93 0.6655
#> 15 15 26 0.6259
#> 16 16 25 0.4586
#> 17 17 85 0.4098
#> 18 18 91 0.2518
#> 19 19 97 0.2148
#> 20 20 95 0.0935
#> 21 21 71 0.0654
#> 22 22 42 0.0568
#> 23 23 31 0.0327
#> 24 24 44 0.0148
#> 25 25 4 0.0083
#> 26 26 15 0.0050
#> 27 27 89 0.0004
#> 28 28 22 0.0001
#> 29 29 58 0.0000
#> 30 30 83 0.0000
#> 31 31 70 0.0000
#> 32 32 64 0.0000
#> 33 33 43 0.0000
#> 34 34 5 0.0000
#> 35 35 6 0.0000
#> 36 36 7 0.0000
#> 37 37 8 0.0000
#> 38 38 10 0.0000
#> 39 39 11 0.0000
#> 40 40 12 0.0000
#> 41 41 13 0.0000
#> 42 42 14 0.0000
#> 43 43 16 0.0000
#> 44 44 19 0.0000
#> 45 45 21 0.0000
#> 46 46 23 0.0000
#> 47 47 24 0.0000
#> 48 48 27 0.0000
#> 49 49 28 0.0000
#> 50 50 29 0.0000head(round(ranking2[["0.85"]], 4), 50)
#> Ranking DMU Probability_predicted betas Probability_target
#> 1 1 1 1.0000 0.0000 0.85
#> 2 2 2 1.0000 0.0000 0.85
#> 3 3 3 1.0000 0.0000 0.85
#> 4 4 17 1.0000 0.0000 0.85
#> 5 5 20 1.0000 0.0000 0.85
#> 6 6 46 1.0000 0.0000 0.85
#> 7 7 9 0.9996 0.0000 0.85
#> 8 8 62 0.9946 0.0000 0.85
#> 9 9 92 0.9768 0.0000 0.85
#> 10 10 18 0.8956 0.0000 0.85
#> 11 11 36 0.8672 0.0000 0.85
#> 12 12 56 0.7489 0.0030 0.85
#> 13 13 75 0.7073 0.0041 0.85
#> 14 14 26 0.6259 0.0094 0.85
#> 15 15 85 0.4098 0.0104 0.85
#> 16 16 25 0.4586 0.0146 0.85
#> 17 17 91 0.2518 0.0164 0.85
#> 18 18 71 0.0654 0.0217 0.85
#> 19 19 31 0.0327 0.0222 0.85
#> 20 20 44 0.0148 0.0247 0.85
#> 21 21 95 0.0935 0.0298 0.85
#> 22 22 42 0.0568 0.0336 0.85
#> 23 23 15 0.0050 0.0348 0.85
#> 24 24 4 0.0083 0.0499 0.85
#> 25 25 70 0.0000 0.0538 0.85
#> 26 26 64 0.0000 0.0666 0.85
#> 27 27 22 0.0001 0.0858 0.85
#> 28 28 58 0.0000 0.0902 0.85
#> 29 29 86 0.0000 0.1415 0.85
#> 30 30 67 0.0000 0.1582 0.85
#> 31 31 48 0.0000 0.1814 0.85
#> 32 32 12 0.0000 0.1819 0.85
#> 33 33 61 0.0000 0.1918 0.85
#> 34 34 6 0.0000 0.2029 0.85
#> 35 35 53 0.0000 0.2084 0.85
#> 36 36 29 0.0000 0.2104 0.85
#> 37 37 68 0.0000 0.2120 0.85
#> 38 38 24 0.0000 0.2245 0.85
#> 39 39 55 0.0000 0.2414 0.85
#> 40 40 37 0.0000 0.2457 0.85
#> 41 41 57 0.0000 0.2460 0.85
#> 42 42 27 0.0000 0.2553 0.85
#> 43 43 60 0.0000 0.2714 0.85
#> 44 44 41 0.0000 0.2731 0.85
#> 45 45 47 0.0000 0.2774 0.85
#> 46 46 35 0.0000 0.2788 0.85
#> 47 47 63 0.0000 0.2870 0.85
#> 48 48 45 0.0000 0.2880 0.85
#> 49 49 14 0.0000 0.3395 0.85
#> 50 50 39 0.0000 0.3464 0.85For each threshold, we identify peer DMUs that serve as reference comparators.
weighted = FALSE treats peers uniformly; set TRUE to
weight importance given relative_importance.
peers <- PEAXAI_peer(
data = data,
x = x,
y = y,
final_model = models[["best_model_fit"]][["nnet"]],
targets = targets,
efficiency_thresholds = efficiency_thresholds,
weighted = FALSE,
relative_importance = relative_importance
)head(peers, 50)
#> DMU 0.75 0.85 0.95
#> 1 1 1 1 1
#> 2 2 2 2 2
#> 3 3 3 3 3
#> 4 4 9 9 9
#> 5 5 9 9 9
#> 6 6 9 9 9
#> 7 7 9 9 9
#> 8 8 3 3 3
#> 9 9 9 9 9
#> 10 10 9 9 9
#> 11 11 9 9 9
#> 12 12 18 18 9
#> 13 13 9 9 9
#> 14 14 9 9 9
#> 15 15 18 18 17
#> 16 16 9 9 9
#> 17 17 17 17 17
#> 18 18 18 18 17
#> 19 19 9 9 9
#> 20 20 20 20 20
#> 21 21 9 9 9
#> 22 22 18 18 17
#> 23 23 9 9 9
#> 24 24 18 18 17
#> 25 25 17 17 17
#> 26 26 17 17 17
#> 27 27 17 17 17
#> 28 28 20 20 20
#> 29 29 18 18 17
#> 30 30 18 18 9
#> 31 31 18 18 17
#> 32 32 18 18 17
#> 33 33 9 9 9
#> 34 34 18 18 17
#> 35 35 18 18 17
#> 36 36 36 36 20
#> 37 37 17 17 17
#> 38 38 18 18 9
#> 39 39 18 18 17
#> 40 40 18 18 46
#> 41 41 17 17 17
#> 42 42 36 36 46
#> 43 43 46 46 46
#> 44 44 46 46 46
#> 45 45 17 17 17
#> 46 46 46 46 46
#> 47 47 17 17 17
#> 48 48 46 46 46
#> 49 49 18 18 17
#> 50 50 46 46 46