Titanic Survival Prediction using Iguanas with Feature Engineering and ONNX file formating.#

This notebook demonstrates a complete end-to-end example of using Iguanas for rule-based classification on the Kaggle Titanic dataset, with advanced feature engineering using the Gators library.

The workflow includes:

Loading and exploring the data
Feature engineering using Gators (new columns, transformations, encoding)
Generating candidate rules using XGBoost
Filtering and selecting high-quality rules
Combining rules using different strategies
Generating predictions for submission

1. Import Libraries#

[1]:

import numpy as np
import polars as pl
from gators.data_cleaning import CastColumns, DropColumns, RenameColumns
from gators.discretizers import CustomDiscretizer
from gators.encoders import RareCategoryEncoder, WOEEncoder
from gators.feature_generation import ConditionFeatures, IsNull, MathFeatures, ScalarMathFeatures
from gators.feature_generation_str import (
    ExtractSubstring,
    SplitExtract,
)
from gators.imputers import NumericImputer, StringImputer
from gators.pipeline import Pipeline
from xgboost import XGBClassifier

from iguanas.metrics import compute_metrics
from iguanas.rule_analysis import generate_rule_performance_report
from iguanas.rule_combination import (
    combine_rules_beam_search,
    combine_rules_cumulative,
    combine_rules_greedy,
)
from iguanas.rule_evaluation import apply_rules
from iguanas.rule_generation import rule_grid_search
from iguanas.rule_selection import filter_correlated_rules

2. Load and Prepare Data#

Load the Titanic training data and separate features from the target variable (Survived).

[2]:

train = pl.read_csv("../../../../../kaggle/titanic/train.csv").drop("PassengerId")
X_train = train.drop("Survived")
y_train = train["Survived"]

3. Feature Engineering with Gators#

Build a comprehensive feature engineering pipeline using the Gators library. This pipeline will:

Create missing value indicators for Age and Cabin
Extract string features (name titles, cabin deck, ticket length)
Calculate family size and fare per person
Create categorical bins for age
Generate feature interactions
Apply Weight of Evidence (WOE) encoding for all categorical variables

Key transformations in this pipeline:

Missing value handling: Create indicators for missing Age/Cabin, then impute
Feature extraction: Extract passenger titles from names, cabin deck letters, ticket lengths
Feature creation: Calculate family size, fare per person, traveling alone indicator
Discretization: Convert continuous Age into categorical bins
Interactions: Create combinations of Pclass, Age, CabinDeck, and Embarked
Encoding: Apply WOE encoding to convert all categorical features to numeric values

[3]:

# Define the feature engineering pipeline
steps = [
    # Create missing value indicators
    ("IsNull", IsNull(subset=["Age", "Cabin"])),
    # String feature engineering
    ("SplitExtractName", SplitExtract(subset=["Name"], by=", ", n=1)),
    ("SplitExtractTitle", SplitExtract(subset=["Name__split_,__1"], by=".", n=0)),
    # Calculate family size (SibSp + Parch + 1)
    (
        "MathFeatures",
        MathFeatures(groups=[["SibSp", "Parch"]], func=["sum"], new_column_names=["Dummy"]),
    ),
    (
        "ScalarMathFeatures",
        ScalarMathFeatures(
            operations=[{"column": "Dummy_sum", "op": "+", "scalar": 1}],
            new_column_names=["FamilySize"],
        ),
    ),
    # Rename for clarity
    (
        "RenameColumns",
        RenameColumns(
            column_mapping={
                "Name__split_,__1__split_._0": "Title",
                # "Cabin__start0_end1": "CabinDeck",
            }
        ),
    ),
    # Handle rare categories (group infrequent values)
    ("RareCategoryEncoder", RareCategoryEncoder(min_count=0.01)),
    # Calculate fare per person
    (
        "MathFeatures2",
        MathFeatures(
            groups=[["Fare", "FamilySize"]], func=["div"], new_column_names=["FarePerPerson"]
        ),
    ),
    # Create 'traveling alone' indicator
    (
        "ConditionFeatures",
        ConditionFeatures(
            conditions=[{"column": "FamilySize", "op": ">", "value": 1}],
            new_column_names=["IsAlone"],
        ),
    ),
    # Drop raw columns no longer needed
    ("DropColumns", DropColumns(subset=["Cabin", "Ticket", "Dummy_sum"])),
    # Impute missing values
    ("NumericImputer", NumericImputer(strategy="mean")),
    ("StringImputer", StringImputer(strategy="constant", value="MISSING")),
    # as_numerics=True outputs float bin indices, required for ONNX export
    ("CustomDiscretizer", CustomDiscretizer(bins={"Age": [0, 12, 18, 35, 60, 100]}, inplace=True, as_numerics=True)),
    # Apply Weight of Evidence encoding (converts all categorical features to numeric)
    ("WOEEncoder", WOEEncoder()),
]

# Build and fit the pipeline
pipe = Pipeline(steps=steps, verbose=True)
X_train_transformed = pipe.fit_transform(X_train, y_train)

print(f"\nOriginal features: {X_train.shape[1]}")
print(f"Engineered features: {X_train_transformed.shape[1]}")

[Pipeline] fit+transform   1/14 · IsNull  |  in: rows=891  cols=10  nulls=866  →  out: rows=891  cols=12  nulls=866  (0.001s)
[Pipeline] fit+transform   2/14 · SplitExtractName  |  in: rows=891  cols=12  nulls=866  →  out: rows=891  cols=12  nulls=866  (0.001s)
[Pipeline] fit+transform   3/14 · SplitExtractTitle  |  in: rows=891  cols=12  nulls=866  →  out: rows=891  cols=12  nulls=866  (0.001s)
[Pipeline] fit+transform   4/14 · MathFeatures  |  in: rows=891  cols=12  nulls=866  →  out: rows=891  cols=13  nulls=866  (0.000s)
[Pipeline] fit+transform   5/14 · ScalarMathFeatures  |  in: rows=891  cols=13  nulls=866  →  out: rows=891  cols=14  nulls=866  (0.000s)
[Pipeline] fit+transform   6/14 · RenameColumns  |  in: rows=891  cols=14  nulls=866  →  out: rows=891  cols=14  nulls=866  (0.000s)
[Pipeline] fit+transform   7/14 · RareCategoryEncoder  |  in: rows=891  cols=14  nulls=866  →  out: rows=891  cols=14  nulls=864  (0.003s)
[Pipeline] fit+transform   8/14 · MathFeatures2  |  in: rows=891  cols=14  nulls=864  →  out: rows=891  cols=15  nulls=864  (0.000s)
[Pipeline] fit+transform   9/14 · ConditionFeatures  |  in: rows=891  cols=15  nulls=864  →  out: rows=891  cols=16  nulls=864  (0.000s)
[Pipeline] fit+transform   10/14 · DropColumns  |  in: rows=891  cols=16  nulls=864  →  out: rows=891  cols=13  nulls=177  (0.000s)
[Pipeline] fit+transform   11/14 · NumericImputer  |  in: rows=891  cols=13  nulls=177  →  out: rows=891  cols=13  nulls=0  (0.001s)
[Pipeline] fit+transform   12/14 · StringImputer  |  in: rows=891  cols=13  nulls=0  →  out: rows=891  cols=13  nulls=0  (0.000s)
[Pipeline] fit+transform   13/14 · CustomDiscretizer  |  in: rows=891  cols=13  nulls=0  →  out: rows=891  cols=13  nulls=0  (0.000s)
[Pipeline] fit+transform   14/14 · WOEEncoder  |  in: rows=891  cols=13  nulls=0  →  out: rows=891  cols=13  nulls=0  (0.004s)

Original features: 10
Engineered features: 13

4. Generate Candidate Rules#

Use XGBoost-based grid search to generate candidate rules from the engineered features. The rule_grid_search_parallel_scales function trains models with different scale_pos_weight values and extracts rules from the decision trees.

[4]:

estimator = XGBClassifier(n_estimators=100, max_depth=4, eval_metric="logloss", random_state=0)
rules = rule_grid_search(
    estimator, X_train_transformed, y_train, scale_pos_weights=np.logspace(0, 3, 50)
)

[5]:

print(f"Number of rules generated: {len(rules)}")

Number of rules generated: 1635

5. Select High-Quality Rules#

Apply the generated rules to the training data, compute performance metrics, and filter based on:

Minimum precision (> 0.15)
Minimum recall (> 0.15)
Maximum correlation between rules (< 0.8)

This ensures we keep only the most useful and diverse rules.

[6]:

R = apply_rules(X_train_transformed, rules.select("rule").to_series().to_list())
M = compute_metrics(R, y_train)
M = M.filter((pl.col("precision") > 0.15) & (pl.col("recall") > 0.15)).sort(
    "accuracy", descending=True
)
importance = dict(zip(M["rule"], M["f0.5"], strict=False))
uncorrelated_rules = filter_correlated_rules(
    R[M["rule"].to_list()], importance=importance, max_corr=0.8
)

[7]:

num_rules = len(uncorrelated_rules)
print(f"Number of selected rules: {num_rules}")

Number of selected rules: 26

6. Combine Rules#

Test different rule combination strategies to find the best performing ruleset.

6.1 Cumulative Combination#

Combines rules cumulatively (rule1 OR rule2 OR … OR ruleN):

[8]:

R_combined = combine_rules_cumulative(
    R[uncorrelated_rules], output_names=[f"combined_rule_{i}" for i in range(1, num_rules + 1)]
)
M_combined = compute_metrics(R_combined, y_train).sort("accuracy", descending=True)
M_combined.head(3)

[8]:

shape: (3, 17)

rule	TP	FP	TN	FN	precision	recall	accuracy	flagged(%)	good_flagged(%)	f0.25	f0.5	f1	f1.5	f2	mcc	num_rules
str	i64	i64	i64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	u32
"combined_rule_4"	252	57	492	90	0.815534	0.736842	0.835017	34.680135	10.382514	0.810443	0.798479	0.774194	0.759388	0.751342	0.646806	1
"combined_rule_5"	252	57	492	90	0.815534	0.736842	0.835017	34.680135	10.382514	0.810443	0.798479	0.774194	0.759388	0.751342	0.646806	1
"combined_rule_6"	252	57	492	90	0.815534	0.736842	0.835017	34.680135	10.382514	0.810443	0.798479	0.774194	0.759388	0.751342	0.646806	1

6.2 Greedy Search#

Uses a greedy algorithm to iteratively select the best rule combination:

[9]:

R_greedy = combine_rules_greedy(R[uncorrelated_rules], y_train, metric="accuracy")
M_greedy = compute_metrics(R_greedy, y_train)
M_greedy

[9]:

shape: (1, 17)

rule	TP	FP	TN	FN	precision	recall	accuracy	flagged(%)	good_flagged(%)	f0.25	f0.5	f1	f1.5	f2	mcc	num_rules
str	i64	i64	i64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	u32
"((X["FamilySize"] < 5.0) & (X[…	270	71	478	72	0.791789	0.789474	0.839506	38.271605	12.932605	0.791652	0.791325	0.79063	0.790185	0.789936	0.660514	4

6.3 Beam Search#

Uses beam search to explore rule combinations up to a maximum number of rules:

[10]:

R_beam = combine_rules_beam_search(R[uncorrelated_rules], y_train, metric="accuracy", max_rules=10)
M_beam = compute_metrics(R_beam, y_train)
M_beam.head(3)

[10]:

shape: (3, 17)

rule	TP	FP	TN	FN	precision	recall	accuracy	flagged(%)	good_flagged(%)	f0.25	f0.5	f1	f1.5	f2	mcc	num_rules
str	i64	i64	i64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	u32
"((X["FamilySize"] < 5.0) & (X[…	270	71	478	72	0.791789	0.789474	0.839506	38.271605	12.932605	0.791652	0.791325	0.79063	0.790185	0.789936	0.660514	4
"((X["FamilySize"] < 5.0) & (X[…	270	71	478	72	0.791789	0.789474	0.839506	38.271605	12.932605	0.791652	0.791325	0.79063	0.790185	0.789936	0.660514	4
"((X["FamilySize"] < 5.0) & (X[…	277	78	471	65	0.780282	0.809942	0.839506	39.842873	14.20765	0.781966	0.786039	0.794835	0.800578	0.803831	0.663424	4

7. Analyze the Best Ruleset#

Generate a detailed report for the best performing ruleset from brute force combination:

[11]:

for r in M_beam["rule"][0].split(" | "):
    print(r)

((X["FamilySize"] < 5.0) & (X["Title"] >= 0.25029))
((X["FarePerPerson_div"] >= 11.5) & (X["Sex"] >= 1.52977))
((X["Cabin__is_null"] < 1.0) & (X["Fare"] >= 7.775) & (X["Fare"] < 151.55) & (X["Age"] < 4.0))
((X["Title"] >= 0.77539) & (X["Pclass"] < 3.0))

[12]:

ruleset = M_beam["rule"][0]
print(f"Selected ruleset: {ruleset}")
report = generate_rule_performance_report(ruleset, X_train_transformed, y_train)
report

Selected ruleset: ((X["FamilySize"] < 5.0) & (X["Title"] >= 0.25029)) | ((X["FarePerPerson_div"] >= 11.5) & (X["Sex"] >= 1.52977)) | ((X["Cabin__is_null"] < 1.0) & (X["Fare"] >= 7.775) & (X["Fare"] < 151.55) & (X["Age"] < 4.0)) | ((X["Title"] >= 0.77539) & (X["Pclass"] < 3.0))

[12]:

shape: (15, 18)

rule_index	rule	TP	FP	TN	FN	precision	recall	accuracy	flagged(%)	good_flagged(%)	f0.25	f0.5	f1	f1.5	f2	mcc	num_rules
str	str	i64	i64	i64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	u32
"0"	"((X["FamilySize"] < 5.0) & (X[…	270	71	478	72	0.791789	0.789474	0.839506	38.271605	12.932605	0.791652	0.791325	0.79063	0.790185	0.789936	0.660514	4
"0.0"	"(X['FamilySize'] < 5.0) & (X['…	239	57	492	103	0.807432	0.69883	0.820426	33.2211	10.382514	0.800118	0.783093	0.749216	0.729	0.718149	0.61435	1
"0.1"	"(X['FarePerPerson_div'] >= 11.…	132	7	542	210	0.94964	0.385965	0.756453	15.600449	1.275046	0.874513	0.734967	0.548857	0.472207	0.437956	0.500197	1
"0.2"	"(X['Cabin__is_null'] < 1.0) & …	72	16	533	270	0.818182	0.210526	0.679012	9.876543	2.91439	0.699429	0.518732	0.334884	0.272886	0.247253	0.295662	1
"0.3"	"(X['Title'] >= 0.77539) & (X['…	166	9	540	176	0.948571	0.48538	0.792368	19.640853	1.639344	0.898154	0.796545	0.642166	0.571202	0.537913	0.574096	1
…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…
"0.2.1"	"(X['Fare'] >= 7.775)"	313	438	111	29	0.416778	0.915205	0.47587	84.287318	79.781421	0.430571	0.467723	0.572736	0.669023	0.738556	0.156873	1
"0.2.2"	"(X['Fare'] < 151.55)"	322	540	9	20	0.37355	0.94152	0.371493	96.74523	98.360656	0.387293	0.424802	0.534884	0.641434	0.721973	-0.115341	1
"0.2.3"	"(X['Age'] < 4.0)"	259	415	134	83	0.384273	0.75731	0.441077	75.645342	75.591985	0.39574	0.426267	0.509843	0.583131	0.634182	0.001575	1
"0.3.0"	"(X['Title'] >= 0.77539)"	249	98	451	93	0.717579	0.72807	0.785634	38.945006	17.850638	0.718188	0.719653	0.722787	0.72481	0.725948	0.548092	1
"0.3.1"	"(X['Pclass'] < 3.0)"	223	177	372	119	0.5575	0.652047	0.667789	44.893378	32.240437	0.562296	0.57415	0.601078	0.619709	0.630656	0.322308	1

[13]:

ruleset

[13]:

'((X["FamilySize"] < 5.0) & (X["Title"] >= 0.25029)) | ((X["FarePerPerson_div"] >= 11.5) & (X["Sex"] >= 1.52977)) | ((X["Cabin__is_null"] < 1.0) & (X["Fare"] >= 7.775) & (X["Fare"] < 151.55) & (X["Age"] < 4.0)) | ((X["Title"] >= 0.77539) & (X["Pclass"] < 3.0))'

8. Generate Predictions on Test Data#

Apply the same preprocessing pipeline to the test data, then use the best ruleset to generate predictions:

[14]:

X_test = pl.read_csv("../../../../../kaggle/titanic/test.csv")
X_test_transformed = pipe.transform(X_test)
y_pred = eval(ruleset.replace("X", "X_test_transformed"))

[Pipeline] transform   1/14 · IsNull  |  in: rows=418  cols=11  nulls=414  →  out: rows=418  cols=13  nulls=414  (0.000s)
[Pipeline] transform   2/14 · SplitExtractName  |  in: rows=418  cols=13  nulls=414  →  out: rows=418  cols=13  nulls=414  (0.001s)
[Pipeline] transform   3/14 · SplitExtractTitle  |  in: rows=418  cols=13  nulls=414  →  out: rows=418  cols=13  nulls=414  (0.001s)
[Pipeline] transform   4/14 · MathFeatures  |  in: rows=418  cols=13  nulls=414  →  out: rows=418  cols=14  nulls=414  (0.000s)
[Pipeline] transform   5/14 · ScalarMathFeatures  |  in: rows=418  cols=14  nulls=414  →  out: rows=418  cols=15  nulls=414  (0.000s)
[Pipeline] transform   6/14 · RenameColumns  |  in: rows=418  cols=15  nulls=414  →  out: rows=418  cols=15  nulls=414  (0.000s)
[Pipeline] transform   7/14 · RareCategoryEncoder  |  in: rows=418  cols=15  nulls=414  →  out: rows=418  cols=15  nulls=414  (0.001s)
[Pipeline] transform   8/14 · MathFeatures2  |  in: rows=418  cols=15  nulls=414  →  out: rows=418  cols=16  nulls=415  (0.000s)
[Pipeline] transform   9/14 · ConditionFeatures  |  in: rows=418  cols=16  nulls=415  →  out: rows=418  cols=17  nulls=415  (0.000s)
[Pipeline] transform   10/14 · DropColumns  |  in: rows=418  cols=17  nulls=415  →  out: rows=418  cols=14  nulls=88  (0.000s)
[Pipeline] transform   11/14 · NumericImputer  |  in: rows=418  cols=14  nulls=88  →  out: rows=418  cols=14  nulls=0  (0.001s)
[Pipeline] transform   12/14 · StringImputer  |  in: rows=418  cols=14  nulls=0  →  out: rows=418  cols=14  nulls=0  (0.000s)
[Pipeline] transform   13/14 · CustomDiscretizer  |  in: rows=418  cols=14  nulls=0  →  out: rows=418  cols=14  nulls=0  (0.000s)
[Pipeline] transform   14/14 · WOEEncoder  |  in: rows=418  cols=14  nulls=0  →  out: rows=418  cols=14  nulls=0  (0.001s)

[15]:

# Create submission file (Kaggle leaderboard score: 0.78)
# Note: +25% better than without feature engineering (0.60)
pl.DataFrame({"PassengerId": X_test["PassengerId"], "Survived": y_pred}).with_columns(
    pl.col("Survived").cast(pl.Int64)
).write_csv("submission_titanic.csv")

[ ]:

9. ONNX Export — Verify End-to-End Predictions Match#

Build two ONNX models and stitch them into a single graph:

Step	Tool	ONNX function
Feature engineering	Gators `Pipeline`	`pipeline_to_onnx`
Rule scoring	Iguanas ruleset	`rules_to_onnx`
End-to-end	combined	`pipeline_to_scoring_onnx`

pipeline_to_scoring_onnx inserts a reshape bridge between the two models: it Unsqueezes each selected preprocessing output column [N] → [N, 1] and Concatenates them into the [N, num_features] matrix that the Iguanas rules model expects as its input X.

[16]:

import onnxruntime as ort
from onnx import TensorProto
from gators.onnx_converters import pipeline_to_onnx, pipeline_to_scoring_onnx

from iguanas.onnx_converter import rules_to_onnx

# --- 1. Iguanas ruleset → ONNX ---
rules_onnx = rules_to_onnx(ruleset)

feature_map = {p.key: p.value for p in rules_onnx.metadata_props}
feature_cols = [feature_map[f"feature_{i}"] for i in range(len(feature_map))]
print(f"Rule features ({len(feature_cols)}): {feature_cols}")

# Ensure pipeline carries the fitted flag (fit_transform sets it from gators ≥ next release)
pipe._is_fitted = True
pipe._input_columns = list(X_train.columns)
pipe._input_dtypes = dict(zip(X_train.columns, X_train.dtypes))

# --- 2. Gators pipeline + Iguanas rules → single end-to-end ONNX model ---
#    pipeline_to_scoring_onnx bridges the preprocessing outputs into the
#    [N, num_features] matrix that rules_onnx expects as its input 'X'.
combined_onnx = pipeline_to_scoring_onnx(pipe, rules_onnx, feature_cols)
print(f"\nCombined ONNX inputs  : {[i.name for i in combined_onnx.graph.input]}")
print(f"Combined ONNX outputs : {[o.name for o in combined_onnx.graph.output]}")

Rule features (8): ['FamilySize', 'Title', 'FarePerPerson_div', 'Sex', 'Cabin__is_null', 'Fare', 'Age', 'Pclass']

Combined ONNX inputs  : ['Pclass__in', 'Name__in', 'Sex__in', 'Age__in', 'SibSp__in', 'Parch__in', 'Ticket__in', 'Fare__in', 'Cabin__in', 'Embarked__in']
Combined ONNX outputs : ['prediction']

[17]:

# --- 3. Build raw input feeds from X_train ---
#    Each pipeline input tensor is named '{col}__in'.
#    Feed dtype must match each input's declared ONNX elem_type exactly.
feeds = {}
for inp in combined_onnx.graph.input:
    col = inp.name.removesuffix("__in")
    s = X_train[col]
    elem_type = inp.type.tensor_type.elem_type
    if elem_type == TensorProto.STRING:
        feeds[inp.name] = s.fill_null("").to_numpy(allow_copy=True)
    elif elem_type == TensorProto.INT64:
        feeds[inp.name] = s.to_numpy(allow_copy=True).astype("int64")
    elif elem_type == TensorProto.DOUBLE:
        feeds[inp.name] = s.to_numpy(allow_copy=True).astype("float64")
    else:
        feeds[inp.name] = s.to_numpy(allow_copy=True).astype("float32")

sess = ort.InferenceSession(combined_onnx.SerializeToString())
onnx_pred = sess.run(None, feeds)[0]  # int64 [N]

# --- 4. Reference predictions via Polars on already-transformed data ---
ref_pred = (
    apply_rules(X_train_transformed, [ruleset])[ruleset]
    .fill_null(False)
    .cast(pl.Int64)
    .to_numpy()
)

match = (onnx_pred == ref_pred).all()
print(f"Predictions match : {match}")
print(f"ONNX   positives  : {int(onnx_pred.sum())} / {len(onnx_pred)}")
print(f"Polars positives  : {int(ref_pred.sum())} / {len(ref_pred)}")
assert match, "ONNX and Polars predictions differ!"

Predictions match : True
ONNX   positives  : 341 / 891
Polars positives  : 341 / 891

[ ]: