# Import file and environment tools.
from pathlib import Path
from importlib.metadata import version
import hashlib
import json
import platform
import joblib

# Your note: 

# Import data tools.
import numpy as np
import pandas as pd

# Your note: 

# Import plotting tools.
import matplotlib.pyplot as plt
import seaborn as sns

# Your note: 

# Import scientific diagnostics.
from scipy import stats

# Your note: 

# Import splitting and validation tools.
from sklearn.model_selection import (
    train_test_split, StratifiedKFold, GridSearchCV, cross_validate
)

# Your note: 

# Import preprocessing tools.
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler, label_binarize

# Your note: 

# Import classifiers.
from sklearn.dummy import DummyClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
from sklearn.tree import DecisionTreeClassifier, plot_tree, export_text

# Your note: 

# Import evaluation tools.
from sklearn.metrics import (
    accuracy_score, balanced_accuracy_score, precision_score, recall_score,
    f1_score, classification_report, confusion_matrix,
    ConfusionMatrixDisplay, roc_curve, roc_auc_score, log_loss
)
from sklearn.inspection import permutation_importance

# Your note: 

# Choose fixed run settings.
RANDOM_SEED = 42
TEST_FRACTION = 0.20
CV_FOLDS = 5
N_JOBS = 1
OUTPUT_DIR = Path('outputs')
OUTPUT_DIR.mkdir(exist_ok=True)

# Your note: 

# Set clear display styles.
pd.set_option('display.max_columns', 12)
pd.set_option('display.max_rows', 30)
pd.set_option('display.precision', 4)
sns.set_theme(style='whitegrid', context='notebook')
plt.rcParams.update({'figure.dpi': 110, 'savefig.dpi': 180})
CLASS_COLORS = ['#667582', '#C7353F', '#247B73', '#16324F', '#AA752C', '#794F91']

# Your note: 

# Record the software environment.
package_versions = {
    name: version(name)
    for name in ['numpy', 'pandas', 'scipy', 'scikit-learn', 'matplotlib', 'seaborn', 'joblib']
}
print('Python:', platform.python_version())
pd.Series(package_versions, name='version')

# Your note: 

# Locate the input file.
DATA_PATH = Path('Stars.csv')
if not DATA_PATH.is_file():
    raise FileNotFoundError('Place Stars.csv beside this notebook, or edit DATA_PATH.')
print('Dataset:', DATA_PATH.name)

# Your note: 

# Fingerprint the input.
data_sha256 = hashlib.sha256(DATA_PATH.read_bytes()).hexdigest()
print('SHA-256:', data_sha256)

# Your note: 

# Load the original data.
raw_df = pd.read_csv(DATA_PATH)

# Your note: 

# Inspect the first five records.
raw_df.head()

# Your note: 

# Count rows and columns.
print(f'{raw_df.shape[0]} rows × {raw_df.shape[1]} columns')

# Your note: 

# Inspect column types.
raw_df.dtypes.rename('dtype')

# Your note: 

# Check missing values once.
raw_df.isna().sum().rename('missing_count')

# Your note: 

# Check duplicate records.
print('Exact duplicate rows:', int(raw_df.duplicated().sum()))

# Your note: 

# Define the class names.
CLASS_NAMES = {
    0: 'Brown Dwarf', 1: 'Red Dwarf', 2: 'White Dwarf',
    3: 'Main Sequence', 4: 'Supergiant', 5: 'Hypergiant'
}
CLASS_IDS = list(CLASS_NAMES)
CLASS_LABELS = list(CLASS_NAMES.values())

# Your note: 

# Define the numeric predictors.
NUMERIC_FEATURES = [
    'Temperature (K)', 'Luminosity (L/Lo)',
    'Radius (R/Ro)', 'Absolute magnitude (Mv)'
]
POSITIVE_FEATURES = NUMERIC_FEATURES[:3]

# Your note: 

# Validate the input schema.
required_columns = set(NUMERIC_FEATURES + ['Star type', 'Star category', 'Star color', 'Spectral Class'])
assert required_columns.issubset(raw_df.columns), 'Required columns are missing.'
assert raw_df['Star type'].notna().all(), 'Missing target labels.'
assert set(raw_df['Star type']) == set(CLASS_IDS), 'Unexpected class codes.'
assert raw_df[NUMERIC_FEATURES].apply(pd.api.types.is_numeric_dtype).all(), 'Non-numeric predictors.'
print('Required columns and class codes are valid.')

# Your note: 

# Verify every target-to-name mapping.
expected_names = raw_df['Star type'].map(CLASS_NAMES)
label_mismatches = raw_df['Star category'].ne(expected_names)
assert not label_mismatches.any(), 'Star category disagrees with Star type.'
print('All numeric class codes match their text labels.')

# Your note: 

# Check numerical validity.
assert np.isfinite(raw_df[NUMERIC_FEATURES].to_numpy()).all(), 'Missing or non-finite measurements.'
assert raw_df[POSITIVE_FEATURES].gt(0).all().all(), 'Non-positive physical measurements.'
print('Predictors are finite; log-transform inputs are positive.')

# Your note: 

# Create a working copy.
df = raw_df.copy(deep=True)

# Your note: 

# Standardise colour spelling.
df['Star color'] = (
    df['Star color'].str.strip().str.lower()
    .str.replace(r'[-\s]+', ' ', regex=True)
)

# Your note: 

# Standardise spectral text.
df['Spectral Class'] = df['Spectral Class'].str.strip().str.upper()

# Your note: 

# Summarise the colour cleaning.
print('Original colour categories:', raw_df['Star color'].nunique())
print('Cleaned colour categories:', df['Star color'].nunique())
df['Star color'].value_counts().rename('records')

# Your note: 

# Summarise the class balance.
class_counts = df['Star type'].value_counts().sort_index()
pd.DataFrame({'class': CLASS_LABELS, 'records': class_counts.to_numpy(),
              'sample_share': class_counts.to_numpy() / len(df)})

# Your note: 

# Separate predictors from labels.
X = df[NUMERIC_FEATURES].copy()
y = df['Star type'].copy()

# Your note: 

# Create the shared stratified split.
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=TEST_FRACTION,
    stratify=y, random_state=RANDOM_SEED
)

# Your note: 

# Verify that records do not cross the split.
assert set(X_train.index).isdisjoint(X_test.index)
assert len(X_train) + len(X_test) == len(df)
assert not X.duplicated().any(), 'Identical feature rows need grouped splitting.'
print(f'Training: {len(X_train)} | Final test: {len(X_test)}')

# Your note: 

# Check the split balance.
split_counts = pd.DataFrame({
    'training': y_train.value_counts().sort_index(),
    'final_test': y_test.value_counts().sort_index()
}).rename(index=CLASS_NAMES)
split_counts

# Your note: 

# Create a training-only exploration table.
train_df = df.loc[X_train.index].copy()

# Your note: 

# Summarise measurement ranges.
train_df[NUMERIC_FEATURES].describe().T

# Your note: 

# Plot training class counts.
fig, ax = plt.subplots(figsize=(9, 4))
ax.bar(CLASS_LABELS, split_counts['training'], color=CLASS_COLORS)
ax.set(title='Training records per dataset class', ylabel='Records', xlabel='Dataset class')
ax.tick_params(axis='x', rotation=25)
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '01_class_balance.png', bbox_inches='tight')
plt.show()

# Your note: 

# Compare typical values within classes.
train_df.groupby('Star category')[NUMERIC_FEATURES].median().reindex(CLASS_LABELS)

# Your note: 

# Plot a Hertzsprung–Russell style diagram.
fig, ax = plt.subplots(figsize=(9, 6))
for class_id, color in zip(CLASS_IDS, CLASS_COLORS):
    group = train_df[train_df['Star type'].eq(class_id)]
    ax.scatter(group['Temperature (K)'], group['Luminosity (L/Lo)'],
               label=CLASS_NAMES[class_id], color=color, s=38, alpha=0.8)
ax.set_xscale('log')
ax.set_yscale('log')
ax.invert_xaxis()
ax.set(title='Training sample: temperature and luminosity', xlabel='Temperature (K; hotter ←)',
       ylabel='Luminosity / solar luminosity (log scale)')
ax.legend(loc='center left', bbox_to_anchor=(1.02, 0.5))
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '02_hr_diagram.png', bbox_inches='tight')
plt.show()

# Your note: 

# Plot temperature against visual magnitude.
fig, ax = plt.subplots(figsize=(9, 6))
sns.scatterplot(data=train_df, x='Temperature (K)', y='Absolute magnitude (Mv)',
                hue='Star category', hue_order=CLASS_LABELS, palette=CLASS_COLORS, ax=ax)
ax.set_xscale('log')
ax.invert_xaxis()
ax.invert_yaxis()
ax.set(title='Training sample: temperature and visual absolute magnitude')
ax.legend(loc='center left', bbox_to_anchor=(1.02, 0.5))
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '03_temperature_magnitude.png', bbox_inches='tight')
plt.show()

# Your note: 

# Inspect feature relationships.
training_correlation = X_train.corr(method='spearman')
fig, ax = plt.subplots(figsize=(8, 6))
sns.heatmap(training_correlation, annot=True, fmt='.2f', cmap='coolwarm',
            vmin=-1, vmax=1, ax=ax)
ax.set_title('Training predictors: Spearman correlation')
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '04_feature_correlation.png', bbox_inches='tight')
plt.show()

# Your note: 

# Create transformed training values for diagnostics.
log_train = X_train.copy()
log_train[POSITIVE_FEATURES] = np.log10(log_train[POSITIVE_FEATURES])

# Your note: 

# Compare transformed distributions.
plot_df = log_train.assign(category=y_train.map(CLASS_NAMES))
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for feature, ax in zip(NUMERIC_FEATURES, axes.flat):
    sns.boxplot(data=plot_df, x='category', y=feature, order=CLASS_LABELS,
                color='#B8C2CC', ax=ax)
    ax.set_xlabel('')
    ax.set_ylabel(('log10 ' if feature in POSITIVE_FEATURES else '') + feature)
    ax.tick_params(axis='x', rotation=35)
fig.suptitle('Training distributions after the chosen transformation', y=1.02)
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '05_transformed_distributions.png', bbox_inches='tight')
plt.show()

# Your note: 

# Check normality within each class.
normality_rows = []
for class_id in CLASS_IDS:
    for feature in NUMERIC_FEATURES:
        values = log_train.loc[y_train.eq(class_id), feature]
        statistic, p_value = stats.shapiro(values)
        normality_rows.append({'class': CLASS_NAMES[class_id], 'feature': feature,
                               'n': len(values), 'W': statistic, 'p_value': p_value})
normality_table = pd.DataFrame(normality_rows)

# Your note: 

# Account for the 24 normality checks.
normality_table['p_bonferroni'] = (normality_table['p_value'] * len(normality_table)).clip(upper=1)
normality_table['evidence_against_normality'] = normality_table['p_bonferroni'].lt(0.05)
normality_table.round(4)

# Your note: 

# Check relationships within classes.
conditional_rows = []
for class_id in CLASS_IDS:
    corr = log_train.loc[y_train.eq(class_id)].corr(method='spearman').abs()
    pairs = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool)).stack()
    pair = pairs.idxmax()
    conditional_rows.append({'class': CLASS_NAMES[class_id],
                             'feature_1': pair[0], 'feature_2': pair[1],
                             'absolute_spearman': pairs.max()})
pd.DataFrame(conditional_rows).round(3)

# Your note: 

# Run an optional physical consistency diagnostic.
implied_luminosity = (
    X_train['Radius (R/Ro)'] ** 2 * (X_train['Temperature (K)'] / 5772.0) ** 4
)
physics_residual = np.log10(X_train['Luminosity (L/Lo)'] / implied_luminosity)
physics_audit = pd.DataFrame({'class': y_train.map(CLASS_NAMES),
                             'log10_observed_over_implied_L': physics_residual})
physics_audit.groupby('class').agg(['median', 'min', 'max']).reindex(CLASS_LABELS).round(3)

# Your note: 

# Define the log transformer.
def make_log_preprocessor():
    return ColumnTransformer(
        [('log10', FunctionTransformer(np.log10, feature_names_out='one-to-one'), POSITIVE_FEATURES)],
        remainder='passthrough', verbose_feature_names_out=False
    )

# Your note: 

# Build the Naive Bayes pipeline.
nb_pipeline = Pipeline([
    ('log_features', make_log_preprocessor()),
    ('scale', StandardScaler()),
    ('classifier', GaussianNB())
])

# Your note: 

# Build the KNN pipeline.
knn_pipeline = Pipeline([
    ('log_features', make_log_preprocessor()),
    ('scale', StandardScaler()),
    ('classifier', KNeighborsClassifier())
])

# Your note: 

# Build the decision tree pipeline.
tree_pipeline = Pipeline([
    ('classifier', DecisionTreeClassifier(random_state=RANDOM_SEED))
])

# Your note: 

# Build the no-skill baseline.
baseline = DummyClassifier(strategy='most_frequent')

# Your note: 

# Show what log transformation and scaling do.
demo_preprocessor = Pipeline([('log_features', make_log_preprocessor()), ('scale', StandardScaler())])
scaled_training = demo_preprocessor.fit_transform(X_train)
pd.DataFrame(scaled_training, columns=NUMERIC_FEATURES).agg(['mean', 'std']).round(3)

# Your note: 

# Define common validation folds.
tuning_cv = StratifiedKFold(n_splits=CV_FOLDS, shuffle=True, random_state=RANDOM_SEED)

# Your note: 

# Define Naive Bayes candidate settings.
nb_grid = {'classifier__var_smoothing': [1e-11, 1e-9, 1e-7, 1e-5, 1e-3]}

# Your note: 

# Define KNN candidate settings.
knn_grid = {
    'classifier__n_neighbors': [1, 3, 5, 7, 9, 11],
    'classifier__weights': ['uniform', 'distance']
}

# Your note: 

# Define tree candidate settings.
tree_grid = {
    'classifier__max_depth': [2, 3, 4, 5, None],
    'classifier__min_samples_leaf': [1, 2, 4]
}

# Your note: 

# Create one search factory.
def make_search(pipeline, grid, cv):
    return GridSearchCV(
        pipeline, grid, cv=cv, scoring={'macro_f1': 'f1_macro', 'accuracy': 'accuracy'},
        refit='macro_f1', return_train_score=True, n_jobs=N_JOBS, error_score='raise'
    )

# Your note: 

# Tune Naive Bayes.
nb_search = make_search(nb_pipeline, nb_grid, tuning_cv)
nb_search.fit(X_train, y_train)
print('Best settings:', nb_search.best_params_)

# Your note: 

# Tune KNN.
knn_search = make_search(knn_pipeline, knn_grid, tuning_cv)
knn_search.fit(X_train, y_train)
print('Best settings:', knn_search.best_params_)

# Your note: 

# Tune the tree.
tree_search = make_search(tree_pipeline, tree_grid, tuning_cv)
tree_search.fit(X_train, y_train)
print('Best settings:', tree_search.best_params_)

# Your note: 

# Collect training search results.
searches = {'Gaussian Naive Bayes': nb_search, 'KNN': knn_search, 'Decision tree': tree_search}
search_rows = []
for name, search in searches.items():
    i = search.best_index_
    search_rows.append({'model': name, 'CV_macro_F1': search.best_score_,
                        'CV_fold_SD': search.cv_results_['std_test_macro_f1'][i],
                        'training_macro_F1': search.cv_results_['mean_train_macro_f1'][i],
                        'best_parameters': str(search.best_params_)})
search_summary = pd.DataFrame(search_rows).set_index('model')
search_summary.round(4)

# Your note: 

# Inspect KNN tuning.
knn_cv_table = pd.DataFrame(knn_search.cv_results_)
fig, ax = plt.subplots(figsize=(8, 4))
for weight, group in knn_cv_table.groupby('param_classifier__weights'):
    ax.errorbar(group['param_classifier__n_neighbors'].astype(int), group['mean_test_macro_f1'],
                yerr=group['std_test_macro_f1'], marker='o', capsize=3, label=weight)
ax.set(title='KNN: training CV macro F1', xlabel='Number of neighbours', ylabel='Macro F1')
ax.legend(title='Weighting')
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '06_knn_tuning.png', bbox_inches='tight')
plt.show()

# Your note: 

# Define inner and outer folds.
inner_cv = StratifiedKFold(n_splits=3, shuffle=True, random_state=RANDOM_SEED)
outer_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_SEED + 1)

# Your note: 

# Evaluate the Naive Bayes tuning procedure.
nb_nested = cross_validate(
    make_search(nb_pipeline, nb_grid, inner_cv), X_train, y_train,
    cv=outer_cv, scoring={'macro_f1': 'f1_macro', 'accuracy': 'accuracy'},
    n_jobs=N_JOBS, error_score='raise'
)

# Your note: 

# Evaluate the KNN tuning procedure.
knn_nested = cross_validate(
    make_search(knn_pipeline, knn_grid, inner_cv), X_train, y_train,
    cv=outer_cv, scoring={'macro_f1': 'f1_macro', 'accuracy': 'accuracy'},
    n_jobs=N_JOBS, error_score='raise'
)

# Your note: 

# Evaluate the tree tuning procedure.
tree_nested = cross_validate(
    make_search(tree_pipeline, tree_grid, inner_cv), X_train, y_train,
    cv=outer_cv, scoring={'macro_f1': 'f1_macro', 'accuracy': 'accuracy'},
    n_jobs=N_JOBS, error_score='raise'
)

# Your note: 

# Compare nested validation results.
nested_results = {'Gaussian Naive Bayes': nb_nested, 'KNN': knn_nested, 'Decision tree': tree_nested}
nested_fold_scores = pd.DataFrame({name: result['test_macro_f1'] for name, result in nested_results.items()})
nested_summary = pd.DataFrame({
    'nested_macro_F1': nested_fold_scores.mean(),
    'outer_fold_SD': nested_fold_scores.std(ddof=1),
    'nested_accuracy': {name: result['test_accuracy'].mean() for name, result in nested_results.items()}
})
nested_summary.index.name = 'model'
nested_summary.sort_values('nested_macro_F1', ascending=False).round(4)

# Your note: 

# Select the family before opening the holdout.
winner_name = nested_summary['nested_macro_F1'].idxmax()
selected_model = searches[winner_name].best_estimator_
print('Selected family:', winner_name)
print('Selected parameters:', searches[winner_name].best_params_)

# Your note: 

# Fit the baseline on training labels.
baseline.fit(X_train, y_train)

# Your note: 

# Collect fitted models.
fitted_models = {'Most-frequent baseline': baseline}
fitted_models.update({name: search.best_estimator_ for name, search in searches.items()})

# Your note: 

# Generate final test predictions.
test_predictions = {name: model.predict(X_test) for name, model in fitted_models.items()}

# Your note: 

# Define classification metrics.
def classification_metrics(actual, predicted):
    return {
        'accuracy': accuracy_score(actual, predicted),
        'balanced_accuracy': balanced_accuracy_score(actual, predicted),
        'macro_precision': precision_score(actual, predicted, labels=CLASS_IDS, average='macro', zero_division=0),
        'macro_recall': recall_score(actual, predicted, labels=CLASS_IDS, average='macro', zero_division=0),
        'macro_F1': f1_score(actual, predicted, labels=CLASS_IDS, average='macro', zero_division=0)
    }

# Your note: 

# Compare the shared holdout metrics.
test_metrics = pd.DataFrame({
    name: classification_metrics(y_test, predicted)
    for name, predicted in test_predictions.items()
}).T
test_metrics.index.name = 'model'
test_metrics['selected_before_test'] = test_metrics.index == winner_name
test_metrics.round(4)

# Your note: 

# Plot model accuracy and macro F1.
fig, ax = plt.subplots(figsize=(9, 5))
test_metrics[['accuracy', 'macro_F1']].plot.bar(ax=ax, color=['#16324F', '#C7353F'])
ax.set(ylim=(0, 1.07), title='Same final holdout: model comparison', ylabel='Score', xlabel='')
ax.tick_params(axis='x', rotation=18)
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '07_model_comparison.png', bbox_inches='tight')
plt.show()

# Your note: 

# Report the selected model result.
selected_predictions = test_predictions[winner_name]
selected_accuracy = test_metrics.loc[winner_name, 'accuracy']
selected_f1 = test_metrics.loc[winner_name, 'macro_F1']
correct_count = int(np.sum(selected_predictions == y_test.to_numpy()))
print(f'{winner_name}: {correct_count}/{len(y_test)} correct ({selected_accuracy:.2%}).')
print(f'Final-test macro F1: {selected_f1:.4f}.')
print(f'Training nested-CV macro F1: {nested_summary.loc[winner_name, "nested_macro_F1"]:.4f}.')
print('These results do not establish accuracy on a new telescope survey.')

# Your note: 

# Define a Wilson accuracy interval.
def wilson_interval(correct, total, confidence=0.95):
    z = stats.norm.ppf(0.5 + confidence / 2)
    proportion = correct / total
    denominator = 1 + z ** 2 / total
    centre = (proportion + z ** 2 / (2 * total)) / denominator
    half_width = z * np.sqrt(proportion * (1 - proportion) / total + z ** 2 / (4 * total ** 2)) / denominator
    return centre - half_width, centre + half_width

# Your note: 

# Quantify holdout uncertainty.
accuracy_low, accuracy_high = wilson_interval(correct_count, len(y_test))
print(f'95% Wilson interval for aggregate accuracy: {accuracy_low:.2%}–{accuracy_high:.2%}')
print(f'One additional test error changes accuracy by {1 / len(y_test):.2%}.')

# Your note: 

# Show all learned-model confusion matrices.
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
for (name, _), ax in zip(searches.items(), axes):
    matrix = confusion_matrix(y_test, test_predictions[name], labels=CLASS_IDS)
    ConfusionMatrixDisplay(matrix, display_labels=CLASS_LABELS).plot(
        ax=ax, cmap='Blues', colorbar=False, xticks_rotation=75
    )
    ax.set_title(name)
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '08_confusion_matrices.png', bbox_inches='tight')
plt.show()

# Your note: 

# Produce the selected model classification report.
report_dict = classification_report(
    y_test, selected_predictions, labels=CLASS_IDS, target_names=CLASS_LABELS,
    output_dict=True, zero_division=0
)
classification_report_table = pd.DataFrame(report_dict).T
classification_report_table.round(4)

# Your note: 

# Build a record-level prediction audit.
prediction_audit = X_test.copy()
prediction_audit.index.name = 'csv_row_index'
prediction_audit['actual_class'] = y_test.map(CLASS_NAMES)
prediction_audit['predicted_class'] = pd.Series(selected_predictions, index=X_test.index).map(CLASS_NAMES)
prediction_audit['correct'] = prediction_audit['actual_class'].eq(prediction_audit['predicted_class'])

# Your note: 

# Inspect the selected model mistakes.
mistakes = prediction_audit.loc[~prediction_audit['correct']]
print('Observed holdout errors:', len(mistakes))
mistakes

# Your note: 

# Compare mistakes from the learned models.
comparison_errors = []
for name in searches:
    for row_index, actual, predicted in zip(X_test.index, y_test, test_predictions[name]):
        if actual != predicted:
            comparison_errors.append({'model': name, 'csv_row_index': row_index,
                                      'actual': CLASS_NAMES[int(actual)],
                                      'predicted': CLASS_NAMES[int(predicted)]})
pd.DataFrame(comparison_errors, columns=['model', 'csv_row_index', 'actual', 'predicted'])

# Your note: 

# Align probability columns explicitly.
raw_probabilities = selected_model.predict_proba(X_test)
probability_order = [list(selected_model.classes_).index(class_id) for class_id in CLASS_IDS]
test_probabilities = raw_probabilities[:, probability_order]
assert np.allclose(test_probabilities.sum(axis=1), 1)

# Your note: 

# Calculate ranking and probability-quality metrics.
test_auc = roc_auc_score(y_test, test_probabilities, labels=CLASS_IDS, multi_class='ovr', average='macro')
test_log_loss = log_loss(y_test, test_probabilities, labels=CLASS_IDS)
print(f'Macro one-vs-rest ROC AUC: {test_auc:.4f}')
print(f'Multiclass log loss (lower is better): {test_log_loss:.4f}')

# Your note: 

# Plot one-vs-rest ROC curves.
binarized_test = label_binarize(y_test, classes=CLASS_IDS)
fig, ax = plt.subplots(figsize=(8, 6))
for i, (class_id, color) in enumerate(zip(CLASS_IDS, CLASS_COLORS)):
    fpr, tpr, _ = roc_curve(binarized_test[:, i], test_probabilities[:, i])
    class_auc = roc_auc_score(binarized_test[:, i], test_probabilities[:, i])
    ax.plot(fpr, tpr, color=color, label=f'{CLASS_NAMES[class_id]} (AUC {class_auc:.3f})')
ax.plot([0, 1], [0, 1], '--', color='#667582', label='Chance ranking')
ax.set(title=f'{winner_name}: final-test one-vs-rest ROC', xlabel='False positive rate', ylabel='True positive rate')
ax.legend(loc='lower right', fontsize=9)
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '09_roc_curves.png', bbox_inches='tight')
plt.show()

# Your note: 

# Attach model scores to the audit table.
prediction_audit['maximum_model_score'] = test_probabilities.max(axis=1)
prediction_audit[['actual_class', 'predicted_class', 'correct', 'maximum_model_score']].head(10)

# Your note: 

# Summarise score bins.
score_bins = pd.cut(prediction_audit['maximum_model_score'],
                    bins=[0, 0.5, 0.8, 0.95, 1.0], include_lowest=True)
score_diagnostic = prediction_audit.groupby(score_bins, observed=True).agg(
    records=('correct', 'size'), mean_score=('maximum_model_score', 'mean'),
    observed_accuracy=('correct', 'mean')
)
score_diagnostic.round(4)

# Your note: 

# Print readable tree rules.
fitted_tree = tree_search.best_estimator_.named_steps['classifier']
print(export_text(fitted_tree, feature_names=NUMERIC_FEATURES, decimals=3))

# Your note: 

# Draw the learned tree.
fig, ax = plt.subplots(figsize=(17, 9))
plot_tree(fitted_tree, feature_names=NUMERIC_FEATURES, class_names=CLASS_LABELS,
          filled=True, rounded=True, impurity=False, proportion=False, ax=ax, fontsize=8)
ax.set_title('Tuned decision tree: rules learned from training records')
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '10_tree_rules.png', bbox_inches='tight')
plt.show()

# Your note: 

# Estimate exploratory feature importance.
importance_result = permutation_importance(
    selected_model, X_train, y_train, scoring='f1_macro',
    n_repeats=20, random_state=RANDOM_SEED, n_jobs=N_JOBS
)
importance_table = pd.DataFrame({
    'feature': NUMERIC_FEATURES, 'mean_F1_drop': importance_result.importances_mean,
    'shuffle_SD': importance_result.importances_std
}).sort_values('mean_F1_drop', ascending=False)
importance_table.round(4)

# Your note: 

# Plot exploratory importance.
fig, ax = plt.subplots(figsize=(9, 4))
ordered_importance = importance_table.sort_values('mean_F1_drop')
ax.barh(ordered_importance['feature'], ordered_importance['mean_F1_drop'],
         xerr=ordered_importance['shuffle_SD'], color='#16324F', capsize=3)
ax.set(title=f'{winner_name}: exploratory training permutation importance',
       xlabel='Drop in training macro F1 after shuffling')
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '11_permutation_importance.png', bbox_inches='tight')
plt.show()

# Your note: 

# Define measurement validation.
def validate_measurements(measurements):
    if not isinstance(measurements, pd.DataFrame):
        raise TypeError('Supply a pandas DataFrame with named measurement columns.')
    if not measurements.columns.is_unique or set(measurements.columns) != set(NUMERIC_FEATURES):
        raise ValueError('Supply exactly the four named predictor columns.')
    if measurements.empty:
        raise ValueError('Supply at least one observation.')
    numeric = measurements.loc[:, NUMERIC_FEATURES].apply(pd.to_numeric, errors='raise')
    if not np.isfinite(numeric.to_numpy()).all():
        raise ValueError('Measurements must be finite and non-missing.')
    if not numeric[POSITIVE_FEATURES].gt(0).all().all():
        raise ValueError('Temperature, luminosity, and radius must be positive.')
    return numeric

# Your note: 

# Record supported marginal ranges.
training_ranges = X_train.agg(['min', 'max'])
training_ranges

# Your note: 

# Define a candidate prediction function.
def suggest_for_review(measurements):
    numeric = validate_measurements(measurements)
    outside = numeric.lt(training_ranges.loc['min']) | numeric.gt(training_ranges.loc['max'])
    supported = ~outside.any(axis=1)
    result = pd.DataFrame(index=numeric.index)
    result['suggested_class'] = 'Abstain'
    result['maximum_model_score'] = np.nan
    result['review_status'] = 'Outside training range; expert review required'
    if supported.any():
        inputs = numeric.loc[supported]
        result.loc[supported, 'suggested_class'] = [CLASS_NAMES[int(c)] for c in selected_model.predict(inputs)]
        result.loc[supported, 'maximum_model_score'] = selected_model.predict_proba(inputs).max(axis=1)
        result.loc[supported, 'review_status'] = 'Candidate only; expert confirmation required'
    return result

# Your note: 

# Create demonstration inputs.
example_measurements = X_test.head(3).copy()
example_measurements

# Your note: 

# Generate candidate suggestions.
example_suggestions = suggest_for_review(example_measurements)
example_suggestions

# Your note: 

# Demonstrate abstention for an unsupported magnitude.
unsupported_example = example_measurements.iloc[[0]].copy()
unsupported_example['Absolute magnitude (Mv)'] = 612.0
suggest_for_review(unsupported_example)

# Your note: 

# Specify workload assumptions.
SCENARIO_RECORDS = 10_000
MANUAL_MINUTES = 3.0
ASSISTED_SAMPLE_RATE = 0.05
INFERENCE_SECONDS = 1.0
ASSUMED_ROUTING_FRACTIONS = [0.0, 0.5, 0.7, 0.9]

# Your note: 

# Calculate hypothetical processing effort.
manual_hours = SCENARIO_RECORDS * MANUAL_MINUTES / 60
scenario_rows = []
for fraction in ASSUMED_ROUTING_FRACTIONS:
    reviewed_fraction = 1 - fraction + fraction * ASSISTED_SAMPLE_RATE
    review_hours = SCENARIO_RECORDS * reviewed_fraction * MANUAL_MINUTES / 60
    processing_hours = SCENARIO_RECORDS * INFERENCE_SECONDS / 3600 if fraction > 0 else 0
    scenario_rows.append({'assumed_routed_share': fraction, 'human_review_hours': review_hours,
                          'serial_processing_hours': processing_hours,
                          'gross_effort_hours': review_hours + processing_hours,
                          'gross_hours_avoided': manual_hours - review_hours - processing_hours})
workload_scenarios = pd.DataFrame(scenario_rows)
workload_scenarios.round(2)

# Your note: 

# Visualise the assumed workload.
fig, ax = plt.subplots(figsize=(8, 4))
ax.bar(workload_scenarios['assumed_routed_share'].map(lambda x: f'{x:.0%}'),
       workload_scenarios['gross_effort_hours'], color='#16324F')
ax.axhline(manual_hours, linestyle='--', color='#C7353F', label='Manual-only assumption')
ax.set(title='Hypothetical workload — not measured business savings',
       xlabel='Assumed share routed through assisted processing', ylabel='Gross effort (hours)')
ax.legend()
fig.tight_layout()
fig.savefig(OUTPUT_DIR / '12_workload_scenario.png', bbox_inches='tight')
plt.show()

# Your note: 

# Export evaluation tables.
test_metrics.to_csv(OUTPUT_DIR / 'model_comparison.csv')
search_summary.to_csv(OUTPUT_DIR / 'training_search_summary.csv')
nested_summary.to_csv(OUTPUT_DIR / 'nested_cv_summary.csv')
nested_fold_scores.to_csv(OUTPUT_DIR / 'nested_cv_fold_scores.csv', index_label='outer_fold')
classification_report_table.to_csv(OUTPUT_DIR / 'selected_classification_report.csv')

# Your note: 

# Export the record-level audit.
prediction_audit.to_csv(OUTPUT_DIR / 'holdout_prediction_audit.csv')
all_model_predictions = pd.DataFrame(test_predictions, index=X_test.index)
all_model_predictions.insert(0, 'actual_class_code', y_test)
all_model_predictions.to_csv(OUTPUT_DIR / 'all_model_holdout_predictions.csv', index_label='csv_row_index')

# Your note: 

# Export diagnostic tables.
normality_table.to_csv(OUTPUT_DIR / 'training_normality_diagnostics.csv', index=False)
physics_audit.to_csv(OUTPUT_DIR / 'training_physics_diagnostic.csv', index_label='csv_row_index')
importance_table.to_csv(OUTPUT_DIR / 'exploratory_training_importance.csv', index=False)
workload_scenarios.to_csv(OUTPUT_DIR / 'hypothetical_workload_scenarios.csv', index=False)
training_ranges.to_csv(OUTPUT_DIR / 'training_input_ranges.csv')

# Your note: 

# Save the complete fitted pipeline.
joblib.dump(selected_model, OUTPUT_DIR / 'selected_stellar_pipeline.joblib')

# Your note: 

# Check that the saved pipeline reproduces predictions.
reloaded_model = joblib.load(OUTPUT_DIR / 'selected_stellar_pipeline.joblib')
assert np.array_equal(reloaded_model.predict(X_test), selected_predictions)
print('Reloaded pipeline reproduces all selected-model holdout predictions.')

# Your note: 

# Build reproducibility metadata.
run_metadata = {
    'project': 'Stellar Classification: From Data to Decisions', 'author': 'Gareth Grant',
    'data_file': DATA_PATH.name, 'data_sha256': data_sha256,
    'source_url_from_original_submission': 'https://www.kaggle.com/datasets/deepu1109/star-dataset',
    'provenance_status': 'Uploaded file verified internally; catalogue lineage and observational metadata unavailable',
    'python_version': platform.python_version(), 'package_versions': package_versions,
    'random_seed': RANDOM_SEED, 'test_fraction': TEST_FRACTION,
    'training_row_indices': sorted(map(int, X_train.index)),
    'test_row_indices': sorted(map(int, X_test.index)),
    'predictors': NUMERIC_FEATURES, 'classes': CLASS_NAMES,
    'selection_rule': 'Highest training nested-CV mean macro F1; fixed family-order tie break',
    'selected_model': winner_name, 'selected_parameters': searches[winner_name].best_params_,
    'final_test_metrics': classification_metrics(y_test, selected_predictions),
    'correct_test_records': correct_count, 'test_records': len(y_test),
    'wilson_accuracy_interval_95': [accuracy_low, accuracy_high],
    'macro_ovr_auc': test_auc, 'multiclass_log_loss': test_log_loss,
    'model_fit_scope': '192 training rows only; final test excluded'
}

# Your note: 

# Save the run metadata.
(OUTPUT_DIR / 'run_metadata.json').write_text(json.dumps(run_metadata, indent=2), encoding='utf-8')
print('Saved run_metadata.json')

# Your note: 

# Generate a result-grounded portfolio summary.
portfolio_summary = (
    '# Stellar classification: a reproducible machine-learning case study\n\n'
    'I rebuilt a six-class classification workflow around 240 labelled stellar records. '
    'The analysis audits the data, visualises physical relationships, and compares Naive Bayes, '
    'KNN, and a decision tree with a no-skill baseline. Learned preprocessing stays inside '
    'training folds; nested cross-validation selects the model family before a shared final test.\n\n'
    f'The selected {winner_name} model classified {correct_count} of {len(y_test)} held-out records '
    f'correctly ({selected_accuracy:.2%}), with macro F1 {selected_f1:.4f}. '
    f'The approximate 95% Wilson accuracy interval is {accuracy_low:.2%} to {accuracy_high:.2%}. '
    'These are results on a small curated dataset and do not establish accuracy on an independent survey.\n\n'
    'The project demonstrates data auditing, reproducible Python, comparative modelling, '
    'uncertainty reporting, explainable rules, and decision-focused interpretation. '
    'Potential extensions include catalogue triage, expert review tools, and research education. '
    'Operational or financial claims would require representative data and a measured pilot.\n'
)
print(portfolio_summary)

# Your note: 

# Save the short website summary.
(OUTPUT_DIR / 'portfolio_summary.md').write_text(portfolio_summary, encoding='utf-8')

# Your note: