"""
Crafte Sessions analysis 2026 09 29 2026 10 29 - the DataLab pipeline as one script

Generated by DataLab on 2026-10-30 14:01. Reproduces the snapshot features, encoding, scaling, segments v27, the ordered comparison,
with the same settings, so the numbers match the app's.

  column definition v3   feature set v6   segments v27   outcome comparison v7

Run it a cell at a time in Spyder (Ctrl+Enter) and look at each table in the
Variable Explorer as you go. Needs pandas, numpy and scikit-learn.
"""

# %% 0. Settings
# Everything the app chose or you chose is here. Change anything and re-run from that point.
import re
import numpy as np
import pandas as pd

DATA_FILE = 'sessions_analysis-2026-09-29-2026-10-29.csv'   # <- the file you imported; set the full path
READ = dict(kind='csv', sep=',', encoding='utf-8', header=True)

ENTITY = 'browserID'        # the same customer across rows
TIME = 'firstSeen'          # when each row happened
SESSION = None       # groups rows into visits
SNAPSHOTS = 'day'     # End of each active day

# The column definition, for the columns used below
COLUMNS = [
    dict(name='sessionDate', role='feature', dtype='datetime', rollup='last', transform='none', missing='leave', use_in_clustering=False, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='lastSeen', role='feature', dtype='datetime', rollup='last', transform='none', missing='leave', use_in_clustering=False, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='minutes', role='feature', dtype='count', rollup='last', transform='log', missing='zero', use_in_clustering=True, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='pages', role='feature', dtype='count', rollup='last', transform='log', missing='zero', use_in_clustering=True, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='source', role='feature', dtype='categorical', rollup='last', transform='one_hot', missing='missing_cat', use_in_clustering=True, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='landingPage', role='feature', dtype='categorical', rollup='last', transform='one_hot_top', missing='missing_cat', use_in_clustering=True, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='device', role='feature', dtype='categorical', rollup='last', transform='one_hot', missing='missing_cat', use_in_clustering=True, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='newVisitor', role='feature', dtype='binary', rollup='ever', transform='none', missing='zero', use_in_clustering=True, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='ordered', role='outcome', dtype='binary', rollup='ever', transform='none', missing='zero', use_in_clustering=False, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
    dict(name='actions', role='feature', dtype='list', rollup='count_each', transform='log', missing='leave', use_in_clustering=True, params={'sep': ',', 'tidy': False, 'merge': '', 'top_n': 15, 'levels': '', 'exclude': ''}),
]

ACTIVITY = ['rows', 'active_days', 'tenure_days']       # activity features fed to the model
EXCLUDE = ['actions[add_checkout]', 'actions[add_featured]', 'actions[add_grid]', 'actions[add_product]', 'actions[add_quickview]', 'actions[co_address]', 'actions[co_card]', 'actions[co_cardform]', 'actions[co_cardload]', 'actions[co_cardpay]', 'actions[co_fail]', 'actions[co_payment]', 'actions[co_place]', 'actions[co_shipping]', 'active_days', 'landingPage=/checkout', 'landingPage=/product/20-traditional-christmas-pudding-round', 'landingPage=/product/30-traditional-christmas-pudding-log-in-cloth', 'landingPage=/product/340-sticky-date-pecan-pudding-800g-log-in-cloth', 'landingPage=/product/363-sticky-date-delite-bites-160g', 'landingPage=/product/366-chocolate-and-pecan-cake-500g', 'landingPage=/product/378-delite-bite-duo', 'landingPage=/product/379-two-todays-weekly-offer', 'landingPage=/product/382-pecan-chocolate-cake-pack', 'landingPage=/product/475-chocolate-and-orange-christmas-pudding-800g-log-in-cloth', 'landingPage=/product/48-chocolate-and-macadamia-christmas-pudding-800g-log-in-cloth', 'source=Chatgpt.com', 'source=Crafte', 'source=Google', 'source=Yahoo', 'source=com.facebook.katana']        # inputs you left out
SCALING = 'standard'
CORR_THRESHOLD = 0.9
PCA = dict(mode='off', n=5, variance=90.0, fit_rows=200000)
SEGMENTS = dict(k_mode='search', k_min=3, k_max=8, k_fixed=8, covariance='diag', n_init=2, min_rows=1, ambiguity=0.6, novelty_pct=1.0, fit_rows=60000)
METHOD = 'bgmm'
GROUP = dict(col='device', mode='even')   # each snapshot compared with others in its group
LABELS = {                                  # your names, merges and 'treat as unknown'
    'S7': dict(label='', merged_into='S5', unknown=False),
    'S8': dict(label='', merged_into='S5', unknown=False),
}
OUTCOME = 'ordered'
TARGETS = ['S2', 'S3', 'S5']     # the target segments
MODEL = dict(kind='forest', split='later', test_share=0.2, before_only=True, add_segment_inputs=True, balanced=False, alert_pct=1.0)

PANDAS_2 = int(pd.__version__.split('.')[0]) >= 2
WEEKDAYS = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
MONTHS = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
TRUE_FALSE = {'true': 1, 'yes': 1, 'y': 1, 't': 1, '1': 1, 'false': 0, 'no': 0, 'n': 0, 'f': 0, '0': 0}


def to_datetime(s):
    # pandas 2 needs format='mixed' for mixed date formats; pandas 1 doesn't know it
    return pd.to_datetime(s, errors='coerce', format='mixed') if PANDAS_2 else pd.to_datetime(s, errors='coerce')


def split_items(s, sep=','):
    # a list value 'add_grid,co_address' -> ['add_grid', 'co_address']; empty -> []
    return [[] if x != x or x is None else [t.strip() for t in str(x).split(sep) if t.strip()] for x in s]


def ranked_counts(v):
    # how often each value appears, most common first; ties in alphabetical order
    vc = v.value_counts()
    if len(vc) < 2:
        return vc
    return vc.iloc[np.lexsort((np.array([str(k) for k in vc.index]), -vc.values))]


def as_number(s):
    # numbers stay numbers; yes/no, true/false and 1/0 text become 1/0; anything else is coerced
    if pd.api.types.is_bool_dtype(s) or pd.api.types.is_numeric_dtype(s):
        return s.astype(float)
    low = s.astype(str).str.strip().str.lower()
    if set(low.dropna().unique()) <= set(TRUE_FALSE) | {'nan'}:
        return low.map(TRUE_FALSE).astype(float)
    return pd.to_numeric(s, errors='coerce')


# %% 1. Load the file
if READ['kind'] == 'excel':
    raw = pd.read_excel(DATA_FILE, header=0 if READ['header'] else None, engine='openpyxl')
else:
    with open(DATA_FILE, encoding=READ['encoding'], errors='replace', newline='') as fh:
        raw = pd.read_csv(fh, sep=READ['sep'], header=0 if READ['header'] else None, low_memory=False)
if not READ['header']:
    raw.columns = ['col_{}'.format(i + 1) for i in range(raw.shape[1])]
raw.columns = [str(c).strip() for c in raw.columns]
print(raw.shape, 'rows x columns')


# %% 2. Clean and order the rows
# Keep only the columns used, drop rows with no entity (or no readable time), then sort each
# entity's rows by time. Everything after this works 'so far' down these sorted rows.
rows = raw[['browserID', 'firstSeen', 'sessionDate', 'lastSeen', 'minutes', 'pages', 'source', 'landingPage', 'device', 'newVisitor', 'ordered', 'actions']].copy()
rows = rows[rows[ENTITY].notna()]
rows['_t'] = to_datetime(rows[TIME])
rows = rows[rows['_t'].notna()]
rows['_order'] = np.arange(len(rows))              # ties in time keep file order
rows = rows.sort_values([ENTITY, '_t', '_order'], kind='mergesort').reset_index(drop=True)
g = pd.Series(pd.factorize(rows[ENTITY].astype(str))[0], index=rows.index)   # entity number for each row

# Which rows end a snapshot
rows['_day'] = rows['_t'].dt.normalize()
snapshot_end = ~rows.duplicated([ENTITY, '_day'], keep='last')       # last row of each active day
entity_last = ~rows.duplicated([ENTITY], keep='last')
print(len(rows), 'rows,', int(snapshot_end.sum()), 'snapshots,', g.nunique(), 'entities')


# %% 3. Running values per entity
# Each value is worked out at every row using only that entity's rows up to it.
rolled, outcome_sofar = {}, {}

# ---- activity, from the ID / time / session columns
rows_so_far = pd.Series(1.0, index=rows.index).groupby(g).cumsum()
first_time = rows.groupby(g)['_t'].transform('min')
tenure_days = (rows['_t'] - first_time).dt.total_seconds() / 86400.0
rolled['rows'] = rows_so_far                        # always kept: decides 'too little history'
rolled['active_days'] = (~pd.DataFrame({'g': g, 'v': rows['_t'].dt.normalize()}).duplicated()).astype(np.int64).groupby(g).cumsum().astype(float)
rolled['tenure_days'] = tenure_days

# ---- sessionDate (feature, datetime): latest value so far
t = to_datetime(rows['sessionDate'])
rolled['sessionDate'] = t.groupby(g).ffill()

# ---- lastSeen (feature, datetime): latest value so far
t = to_datetime(rows['lastSeen'])
rolled['lastSeen'] = t.groupby(g).ffill()

# ---- minutes (feature, count): latest value so far; missing values: zero
x = as_number(rows['minutes'])
x = x.fillna(0)
x = x.astype(float)
rolled['minutes'] = x.groupby(g).ffill()

# ---- pages (feature, count): latest value so far; missing values: zero
x = as_number(rows['pages'])
x = x.fillna(0)
x = x.astype(float)
rolled['pages'] = x.groupby(g).ffill()

# ---- source (feature, categorical): latest value so far; missing values: missing_cat
v = rows['source'].astype(object).where(rows['source'].notna())
v = v.fillna('missing')
v = v.where(v.isna(), v.astype(str))
rolled['source'] = v.groupby(g).ffill()

# ---- landingPage (feature, categorical): latest value so far; missing values: missing_cat
v = rows['landingPage'].astype(object).where(rows['landingPage'].notna())
v = v.fillna('missing')
v = v.where(v.isna(), v.astype(str))
rolled['landingPage'] = v.groupby(g).ffill()

# ---- device (feature, categorical): latest value so far; missing values: missing_cat
v = rows['device'].astype(object).where(rows['device'].notna())
v = v.fillna('missing')
v = v.where(v.isna(), v.astype(str))
rolled['device'] = v.groupby(g).ffill()

# ---- newVisitor (feature, binary): ever true so far; missing values: zero
x = as_number(rows['newVisitor'])
x = x.fillna(0)
x = x.astype(float)
rolled['newVisitor'] = x.groupby(g).cummax().groupby(g).ffill()

# ---- ordered (outcome, binary): ever true so far; missing values: zero
x = as_number(rows['ordered'])
x = x.fillna(0)
x = x.astype(float)
outcome_sofar['ordered'] = x.groupby(g).cummax().groupby(g).ffill()

# ---- actions (feature, list): count of each value so far
items = split_items(rows['actions'], sep=',')        # each row's items, e.g. 'a,b' -> ['a', 'b']
counts = ranked_counts(pd.Series([x for r in items for x in r], dtype=object))
values = [str(k) for k in counts.index[:len(counts) if len(counts) <= 60 else 15]]
for val in values:                                   # one running count per item
    rolled['actions[' + val + ']'] = pd.Series([float(r.count(val)) for r in items], index=rows.index).groupby(g).cumsum()
if len(counts) > len(values):                        # everything else, together
    rolled['actions[other]'] = pd.Series([float(sum(1 for x in r if x not in values)) for r in items], index=rows.index).groupby(g).cumsum()


# %% 4. The snapshot table: one row per snapshot
snaps = pd.DataFrame({'entity': rows.loc[snapshot_end, ENTITY].astype(str).values,
                      'snapshot_time': rows.loc[snapshot_end, '_t'].values})
snaps['snapshot_no'] = snaps.groupby('entity').cumcount() + 1
for name, series in rolled.items():
    snaps[name] = series[snapshot_end].values
# outcomes: 'so far' at the snapshot, and 'eventually' at the entity's last row
entity_of_snapshot = g[snapshot_end].values
for name, series in outcome_sofar.items():
    final = pd.Series(series[entity_last].values, index=g[entity_last].values)
    snaps['y_sofar:' + name] = series[snapshot_end].values
    snaps['y_final:' + name] = final.reindex(entity_of_snapshot).values
snaps.head()


# %% 5. Encode the inputs the model sees
inputs = {}
inputs['rows'] = np.log1p(snaps['rows'].astype(float))      # activity: log
inputs['active_days'] = np.log1p(snaps['active_days'].astype(float))      # activity: log
inputs['tenure_days'] = np.log1p(snaps['tenure_days'].astype(float))      # activity: log

# ---- minutes: log (compresses long tails)
for col in ['minutes']:
    x = snaps[col].astype(float)
    x = np.sign(x) * np.log1p(x.abs()) if (x.dropna() < 0).any() else np.log1p(x)
    inputs[col] = x

# ---- pages: log (compresses long tails)
for col in ['pages']:
    x = snaps[col].astype(float)
    x = np.sign(x) * np.log1p(x.abs()) if (x.dropna() < 0).any() else np.log1p(x)
    inputs[col] = x

# ---- source: one 0/1 input per value
v = snaps['source'].astype(object).where(snaps['source'].notna(), 'missing').astype(str)
counts = ranked_counts(v)
values = [str(k) for k in counts.index[:60]]
for val in values:
    inputs['source' + '=' + val] = (v == val).astype(float)
if len(counts) > len(values):
    inputs['source=other'] = (~v.isin(values)).astype(float)

# ---- landingPage: one 0/1 input for each common value, plus 'other'
v = snaps['landingPage'].astype(object).where(snaps['landingPage'].notna(), 'missing').astype(str)
counts = ranked_counts(v)
values = [str(k) for k in counts.index[:15]]
for val in values:
    inputs['landingPage' + '=' + val] = (v == val).astype(float)
if len(counts) > len(values):
    inputs['landingPage=other'] = (~v.isin(values)).astype(float)

# ---- device: one 0/1 input per value
v = snaps['device'].astype(object).where(snaps['device'].notna(), 'missing').astype(str)
counts = ranked_counts(v)
values = [str(k) for k in counts.index[:60]]
for val in values:
    inputs['device' + '=' + val] = (v == val).astype(float)
if len(counts) > len(values):
    inputs['device=other'] = (~v.isin(values)).astype(float)

# ---- newVisitor: used as is
for col in ['newVisitor']:
    x = snaps[col].astype(float)
    inputs[col] = x

# ---- actions: log (compresses long tails)
for col in [k for k in snaps.columns if k.startswith('actions[')]:
    x = snaps[col].astype(float)
    x = np.sign(x) * np.log1p(x.abs()) if (x.dropna() < 0).any() else np.log1p(x)
    inputs[col] = x

X = pd.DataFrame(inputs, index=snaps.index)

# fill anything still empty (e.g. a mean before any value was seen) with that input's median
for col in X.columns:
    if X[col].isna().any():
        X[col] = X[col].fillna(X[col].median() if X[col].notna().any() else 0.0)

# inputs that never vary can't separate anyone; EXCLUDE are the ones you left out
sd = X.std(ddof=0)
constant = [c for c in X.columns if not np.isfinite(sd[c]) or sd[c] == 0]
used = [c for c in X.columns if c not in constant and c not in EXCLUDE]
print('inputs:', len(used), 'used of', X.shape[1], '| never vary:', constant)
Xu = X[used]


# %% 6. Scale, check correlations, optional PCA
if SCALING == 'standard':
    center, scale = Xu.mean(), Xu.std(ddof=0)
elif SCALING == 'robust':
    center = Xu.median()
    scale = Xu.quantile(0.75) - Xu.quantile(0.25)
    scale = scale.where(scale > 0, Xu.std(ddof=0))
elif SCALING == 'minmax':
    center, scale = Xu.min(), Xu.max() - Xu.min()
else:
    center, scale = pd.Series(0.0, index=used), pd.Series(1.0, index=used)
Z = (Xu - center) / scale.where(scale > 0, 1.0)      # the scaled inputs

corr = Z.corr()
pairs = (corr.where(np.triu(np.ones(corr.shape, bool), 1)).stack().rename('r').reset_index()
         .assign(abs_r=lambda t: t['r'].abs()).sort_values('abs_r', ascending=False))
print(pairs[pairs['abs_r'] >= CORR_THRESHOLD].head(20).to_string(index=False))

from sklearn.decomposition import PCA as _PCA
full_pca = _PCA(n_components=int(min(Z.shape[1], len(Z), 50)), random_state=0).fit(
    Z.sample(PCA['fit_rows'], random_state=0).values if len(Z) > PCA['fit_rows'] else Z.values)
explained = pd.Series(full_pca.explained_variance_ratio_.cumsum(), index=range(1, full_pca.n_components_ + 1))
print('components needed for 80% / 90%:', int((explained < 0.8).sum() + 1), '/', int((explained < 0.9).sum() + 1))

if PCA['mode'] == 'off':
    M = Z.reset_index(drop=True)                      # the segmentation uses the inputs directly
else:
    k_all = full_pca.n_components_
    if PCA['mode'] == 'variance':
        target = max(1.0, min(100.0, float(PCA['variance']))) / 100
        k = int(np.searchsorted(explained.values, target) + 1) if explained.values[-1] >= target else k_all
    else:
        k = max(1, min(int(PCA['n']), k_all))
    sample = Z.sample(PCA['fit_rows'], random_state=0) if len(Z) > PCA['fit_rows'] else Z
    pca = _PCA(n_components=k, random_state=0).fit(sample.values)
    M = pd.DataFrame(pca.transform(Z.values), columns=['PC{}'.format(i + 1) for i in range(k)])
    loadings = pd.DataFrame(pca.components_, columns=used, index=M.columns).T   # what drives each component
print('model matrix:', M.shape)


# %% 7. Segments: Gaussian mixture, best number of segments by BIC
from sklearn.mixture import GaussianMixture
# group handling for device: its own inputs are taken out of the segmentation
if GROUP['col'] in snaps.columns:
    groups = snaps[GROUP['col']].astype(object).where(snaps[GROUP['col']].notna(), '(none)').astype(str).values
else:                                               # counted per value: the most frequent so far
    cnt = [c for c in snaps.columns if c.startswith(GROUP['col'] + '[') and c.endswith(']')]
    sub = snaps[cnt].apply(pd.to_numeric, errors='coerce').fillna(0)
    best = np.array([c[len(GROUP['col']) + 1:-1] for c in cnt], dtype=object)[sub.values.argmax(axis=1)]
    groups = np.where(sub.values.max(axis=1) > 0, best, '(none)').astype(str)
own = [c for c in M.columns if c == GROUP['col'] or c.startswith(GROUP['col'] + '=') or c.startswith(GROUP['col'] + '[')]
M = M.drop(columns=own)

# even out: standardise each input within each group, then put it back on the overall scale,
# so a desktop snapshot is 'high' only if it's high for a desktop (groups under 20 left as they are)
def even_out(df, groups, min_n=20):
    out = df.copy().astype(float)
    mu, sd = out.mean(), out.std(ddof=0).replace(0, 1.0)
    g = pd.Series(groups, index=out.index)
    for val, idx in g.groupby(g).groups.items():
        if len(idx) < min_n:
            continue
        part = out.loc[idx]
        gsd = part.std(ddof=0)
        gsd = gsd.where(gsd > 0, sd)
        out.loc[idx] = (part - part.mean()) / gsd * sd + mu
    return out
M = even_out(M, groups)
own_z = [c for c in Z.columns if c == GROUP['col'] or c.startswith(GROUP['col'] + '=') or c.startswith(GROUP['col'] + '[')]
Z = pd.concat([even_out(Z.drop(columns=own_z).reset_index(drop=True), groups), Z[own_z].reset_index(drop=True)], axis=1)[list(Z.columns)]

# the clustering method: Bayesian Gaussian mixture
def fit_method(X, e, S):
    # Bayesian Gaussian mixture: starts with up to k_max segments and switches off the ones it doesn't need
    from sklearn.mixture import BayesianGaussianMixture
    Xe = X[e]
    rng = np.random.RandomState(0)
    fit_idx = rng.choice(len(Xe), S['fit_rows'], replace=False) if len(Xe) > S['fit_rows'] else np.arange(len(Xe))
    Xf = Xe[fit_idx]
    kmax = max(2, min(S['k_max'] if S['k_mode'] == 'search' else S['k_fixed'], len(Xf) - 1))
    bg = BayesianGaussianMixture(n_components=kmax, covariance_type=S['covariance'], n_init=S['n_init'],
                                 weight_concentration_prior_type='dirichlet_process', weight_concentration_prior=1.0 / kmax,
                                 reg_covar=1e-3, max_iter=500, random_state=0).fit(Xf)
    full, ll = np.zeros((len(X), kmax)), np.full(len(X), np.nan)
    full[e], ll[e] = bg.predict_proba(X[e]), bg.score_samples(X[e])
    used = [j for j in range(kmax) if (full[e].argmax(axis=1) == j).any()]    # the segments it kept
    post = full[:, used]
    tot = post.sum(axis=1, keepdims=True)
    post = np.where(tot > 0, post / np.where(tot > 0, tot, 1), 0.0)
    print(len(used), 'segments kept of', kmax)
    return dict(post=post, ll=ll, cut=np.percentile(bg.score_samples(Xf), S['novelty_pct']), noise=None)


Xm = M.values.astype(float)
eligible = (snaps['rows'] >= SEGMENTS['min_rows']).values      # enough history to say anything
r = fit_method(Xm, eligible, SEGMENTS)
post, loglik = r['post'], r['ll']
k = post.shape[1]
raw_best = post.argmax(axis=1)
order = np.argsort(-np.bincount(raw_best[eligible], minlength=k), kind='mergesort')   # S1 = the largest segment; ties keep their order
rename = {int(old): 'S{}'.format(new + 1) for new, old in enumerate(order)}
nearest = np.array([rename[int(b)] for b in raw_best], dtype=object)
confidence = post.max(axis=1)
state = nearest.copy()
state[confidence < SEGMENTS['ambiguity']] = 'U-ambiguous'      # between segments
state[loglik < r['cut']] = 'U-novel'                           # fits no segment
if r['noise'] is not None:
    state[r['noise']] = 'U-novel'                              # density methods: their own 'noise'
state[~eligible] = 'U-insufficient'                            # too little history
assign = pd.DataFrame({'state': state, 'confidence': np.where(eligible, confidence, np.nan), 'loglik': loglik})
for old, new in rename.items():
    assign['p_' + new] = post[:, old]
print(assign['state'].value_counts().sort_index())

# your labels and merges: map each segment to the group it now belongs to
group_of = {s: s for s in set(state)}
for s, lab in LABELS.items():
    if lab.get('unknown') and not s.startswith('U-'):
        group_of[s] = 'U-manual'
    elif lab.get('merged_into') and lab['merged_into'] in group_of and lab['merged_into'] != s:
        group_of[s] = lab['merged_into']
for _ in range(5):                                  # follow chains such as S5 -> S3 -> S2
    group_of = {s: group_of.get(t, t) for s, t in group_of.items()}
assign['group'] = assign['state'].map(group_of)


# %% 8. Describe the segments
# average of each scaled input per segment, in standard deviations from everyone
Zs = (Z - Z.mean()) / Z.std(ddof=0).replace(0, 1)
profile = Zs.groupby(assign['group'].values).mean().T
print(profile.round(1).to_string())
for name in [c for c in snaps.columns if c.startswith('y_final:')]:
    rate = (pd.to_numeric(snaps[name], errors='coerce').fillna(0) > 0).groupby(assign['group'].values).mean()
    print('\n' + name[8:] + ' rate by segment (overall {:.2%}):'.format(rate.mean() if False else
          (pd.to_numeric(snaps[name], errors='coerce').fillna(0) > 0).mean()))
    print(rate.sort_values(ascending=False).round(4).to_string())

# movement: where each entity's next snapshot lands (rows: from, columns: to)
moves = pd.DataFrame({'e': snaps['entity'], 'n': snaps['snapshot_no'], 's': assign['group']}).sort_values(['e', 'n'], kind='mergesort')
moves['next'] = moves.groupby('e')['s'].shift(-1)
transitions = pd.crosstab(moves['s'], moves['next'], normalize='index').round(2)

# optional heatmap of the profile (uncomment)
# import matplotlib.pyplot as plt
# plt.figure(figsize=(1 + 0.6 * profile.shape[1], 0.3 * profile.shape[0] + 1))
# plt.imshow(profile.clip(-3, 3), cmap='RdBu_r', vmin=-3, vmax=3, aspect='auto')
# plt.xticks(range(profile.shape[1]), profile.columns); plt.yticks(range(profile.shape[0]), profile.index)
# plt.colorbar(label='sd from everyone'); plt.tight_layout(); plt.show()


# %% 9. Outcome: training and test sets
y = pd.to_numeric(snaps['y_final:' + OUTCOME], errors='coerce').fillna(0).values > 0      # eventually
y_sofar = pd.to_numeric(snaps['y_sofar:' + OUTCOME], errors='coerce').fillna(0).values > 0  # already, at the snapshot
group = assign['group'].values

# Is the outcome an event, or a label written on every row from the start? For a label
# there is no 'before', so dropping snapshots after it would drop them all.
first_snap = snaps['snapshot_no'].values == 1
from_start = y_sofar[y & first_snap].mean() if (y & first_snap).any() else 0.0
label_style = from_start >= 0.9
keep = ~y_sofar if (MODEL['before_only'] and not label_style) else np.ones(len(y), bool)
print('label written on every row' if label_style else 'event', '| snapshots kept:', int(keep.sum()))

Xo = Z.reset_index(drop=True).copy()                # the scaled inputs (before any PCA)
if MODEL['add_segment_inputs']:
    for col in [c for c in assign.columns if c.startswith('p_')]:
        Xo['segment ' + col[2:] + ' probability'] = assign[col].fillna(0).values
    Xo['fit to segments'] = assign['loglik'].fillna(assign['loglik'].min()).values
    Xo['segment confidence'] = assign['confidence'].fillna(0).values

# split by entity, so nobody is in both training and test
ent = snaps['entity'].astype(str).values
first_seen = pd.Series(pd.to_datetime(snaps['snapshot_time']).values, index=ent).groupby(level=0).min()
if MODEL['split'] == 'random':
    test_entities = set(first_seen.index[np.random.RandomState(0).rand(len(first_seen)) < MODEL['test_share']])
else:
    cut = first_seen.quantile(1 - MODEL['test_share'])          # later customers
    test_entities = set(first_seen.index[first_seen > cut])
test = np.array([e in test_entities for e in ent])
in_target = np.isin(group, TARGETS)
tr, te = keep & ~test, keep & test
print('training rate: everyone {:.2%} -> target segments {:.2%}; targets hold {:.0%} of the outcome'.format(
    y[tr].mean(), y[tr & in_target].mean(), y[tr & in_target].sum() / max(y[tr].sum(), 1)))


# %% 10. Three approaches, scored on the same test snapshots
from sklearn.ensemble import RandomForestClassifier
make_model = lambda: RandomForestClassifier(n_estimators=200, min_samples_leaf=3, random_state=0, n_jobs=1)
from sklearn.metrics import average_precision_score, roc_auc_score
from sklearn.utils.class_weight import compute_sample_weight
Xv = Xo.values
weights = (lambda yy: compute_sample_weight('balanced', yy)) if MODEL['balanced'] else (lambda yy: None)

# 1) segments only: each snapshot scores its segment's outcome rate in training
rates = pd.Series(y[tr]).groupby(group[tr]).mean()
seg_score = np.array([rates.get(s, y[tr].mean()) for s in group])

# 2) reduced: a model trained only on the target segments; everyone else ranks below them
m_reduced = make_model().fit(Xv[tr & in_target], y[tr & in_target], sample_weight=weights(y[tr & in_target]))
red_score = np.where(in_target, 1.0 + m_reduced.predict_proba(Xv)[:, 1], seg_score)

# 3) full: the same model trained on everything
m_full = make_model().fit(Xv[tr], y[tr], sample_weight=weights(y[tr]))
full_score = m_full.predict_proba(Xv)[:, 1]

alert_n = max(1, int(round(te.sum() * MODEL['alert_pct'] / 100.0)))
results = []
for name, score in (('segments only', seg_score), ('model on target segments', red_score), ('model on everything', full_score)):
    yt, st = y[te], score[te]
    top = yt[np.argsort(-st, kind='mergesort')[:alert_n]]
    results.append(dict(approach=name, pr_auc=average_precision_score(yt, st), roc_auc=roc_auc_score(yt, st),
                        caught=int(top.sum()), of=int(yt.sum()), precision=top.mean(), recall=top.sum() / yt.sum()))
results = pd.DataFrame(results)
print('top {} test snapshots flagged'.format(alert_n))
print(results.round(3).to_string(index=False))


# %% 11. Is the difference more than chance?
# Re-score both models on 200 resamples of the test entities (whole entities, since one entity's
# snapshots aren't independent). 'Better' only if the reduced model is ahead in 95% of them.
yt, a, b, et = y[te], red_score[te], full_score[te], ent[te]
codes, uniq = pd.factorize(et)
rows_of = [np.where(codes == i)[0] for i in range(len(uniq))]
rng = np.random.RandomState(0)
diffs = []
for _ in range(200):
    idx = np.concatenate([rows_of[i] for i in rng.randint(0, len(rows_of), len(rows_of))])
    if 0 < yt[idx].sum() < len(idx):
        diffs.append(average_precision_score(yt[idx], a[idx]) - average_precision_score(yt[idx], b[idx]))
lo, hi = np.percentile(diffs, [2.5, 97.5])
print('PR-AUC difference (target segments minus everything): {:+.3f}, 95% range {:+.3f} to {:+.3f}'.format(
    average_precision_score(yt, a) - average_precision_score(yt, b), lo, hi))
print('verdict:', 'better' if lo > 0 else ('worse' if hi < 0 else 'about the same - the gap could be chance'))


# %% 12. What the target-segment model uses (slower; optional)
from sklearn.inspection import permutation_importance
mask = te & in_target if (y[te & in_target].sum() >= 5 and (~y[te & in_target]).sum() >= 5) else tr & in_target
imp = permutation_importance(m_reduced, Xv[mask], y[mask], scoring='average_precision', n_repeats=3, random_state=0, n_jobs=1)
importance = pd.Series(imp.importances_mean, index=Xo.columns).sort_values(ascending=False)
print(importance.head(15).round(4).to_string())