Updated fanny algorithm, fuzzy clustering, and added property based clustering

Author: yuqi-liang-qiqi

.gitignore

@@ -59,6 +59,7 @@ developer/WeightedCluster-master/
 developer/seqHMM-main/
 developer/Py_FS-main/
 developer/eBoruta-master/
+developer/cluster-master/
 
 # ignore RData and Rhistory
 .RData

Tutorials/cluster_analysis/test_sequences_to_variables.py

@@ -131,7 +131,7 @@ def test_phase2(diss_small, kmed, k):
 
 def test_phase3(diss_small, k):
     """Phase 3: fanny_membership and soft_classification_variables."""
-    U, medoids = fanny_membership(diss_small, k=k, m=1.4, random_state=42)
+    U, medoids = fanny_membership(diss_small, k=k, m=1.4)
     n_use = diss_small.shape[0]
     assert U.shape == (n_use, k)
     np.testing.assert_allclose(U.sum(axis=1), 1.0, rtol=1e-5)

Tutorials/cluster_analysis/test_studer_2018_workflow.py

@@ -0,0 +1,177 @@
+#!/usr/bin/env python3
+"""
+Studer (2018) fuzzy + property-based clustering workflow (biofam).
+
+Mirrors the WeightedCluster short R tutorial:
+  seqdef -> seqdist(LCS) -> fanny -> summary(membership)
+  -> seqpropclust(state, duration) -> as.clustrange
+
+Run from repo root:
+    python3 -u Tutorials/cluster_analysis/test_studer_2018_workflow.py
+    python3 -u Tutorials/cluster_analysis/test_studer_2018_workflow.py --full
+
+Default uses n=400 sequences (fast sanity check). --full uses all 2000 rows;
+FANNY alone may take 30-60+ minutes in pure Python (R cluster::fanny ~35 s).
+"""
+from __future__ import annotations
+
+import argparse
+import sys
+import time
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+_script_dir = Path(__file__).resolve().parent
+_repo_root = _script_dir.parents[1]
+if str(_repo_root) not in sys.path:
+    sys.path.insert(0, str(_repo_root))
+
+from sequenzo import SequenceData, get_distance_matrix, load_dataset
+from sequenzo.clustering import (
+    get_fuzzy_clusters,
+    membership_summary,
+    print_property_tree,
+    property_based_clustering,
+    property_clustering_quality,
+)
+
+# R tutorial reference targets (WeightedCluster on biofam, set.seed(1))
+R_MEMBERSHIP_SUMMARY_MEAN = {
+    "V1": 0.205773,
+    "V2": 0.211530,
+    "V3": 0.210022,
+    "V4": 0.178295,
+    "V5": 0.194380,
+}
+R_PROPERTY_GLOBAL_R2 = 0.48761
+R_CLUSTER_QUALITY = pd.DataFrame(
+    {
+        "PBC": [0.50, 0.54, 0.54, 0.58],
+        "HG": [0.61, 0.69, 0.73, 0.80],
+        "R2": [0.21, 0.34, 0.42, 0.49],
+    },
+    index=["cluster2", "cluster3", "cluster4", "cluster5"],
+)
+
+STATE_LABELS = [
+    "Parent",
+    "Left",
+    "Married",
+    "Left/Married",
+    "Child",
+    "Left/Child",
+    "Left/Married/Child",
+    "Divorced",
+]
+TIME_COLS = [str(age) for age in range(15, 31)]
+
+
+def load_biofam_sequences(n_rows: int | None) -> SequenceData:
+    df = load_dataset("biofam").reset_index(drop=True)
+    if n_rows is not None:
+        df = df.head(n_rows)
+    df["id"] = np.arange(len(df))
+    return SequenceData(
+        df,
+        time=TIME_COLS,
+        states=list(range(8)),
+        labels=STATE_LABELS,
+        id_col="id",
+    )
+
+
+def _compare_table(label: str, got: pd.DataFrame, ref: pd.DataFrame, cols: list[str]) -> None:
+    print(f"\n=== {label} (Sequenzo vs R reference) ===")
+    for col in cols:
+        if col not in got.columns or col not in ref.columns:
+            continue
+        g = got[col].to_numpy(dtype=float)
+        r = ref[col].to_numpy(dtype=float)
+        diff = np.abs(g - r)
+        print(f"  {col}: max |diff| = {diff.max():.4f}, mean |diff| = {diff.mean():.4f}")
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Studer 2018 biofam workflow test")
+    parser.add_argument(
+        "--full",
+        action="store_true",
+        help="Use all 2000 biofam sequences (slow FANNY)",
+    )
+    parser.add_argument(
+        "--n",
+        type=int,
+        default=400,
+        help="Subsample size when not using --full (default: 400)",
+    )
+    args = parser.parse_args()
+    n_rows = None if args.full else args.n
+
+    print(f"[>] Loading biofam ({'full' if args.full else f'n={n_rows}'})...")
+    t0 = time.perf_counter()
+    seqdata = load_biofam_sequences(n_rows)
+    print(f"    {seqdata.seqdata.shape[0]} sequences, {len(TIME_COLS)} time points")
+
+    print("[>] LCS distance matrix...")
+    t1 = time.perf_counter()
+    diss = get_distance_matrix(seqdata, method="LCS")
+    diss = np.asarray(diss.values if hasattr(diss, "values") else diss, dtype=float)
+    print(f"    done in {time.perf_counter() - t1:.1f}s; mean diss = {diss.mean():.4f}")
+
+    print("[>] FANNY fuzzy clustering (k=5, memb.exp=1.5)...")
+    t2 = time.perf_counter()
+    fclust = get_fuzzy_clusters(diss, n_clusters=5, memb_exp=1.5, method="fanny")
+    print(
+        f"    done in {time.perf_counter() - t2:.1f}s; "
+        f"converged={fclust.converged}, iterations={fclust.iterations}"
+    )
+
+    summary = membership_summary(fclust.membership)
+    print("\n=== membership summary (Mean row) ===")
+    print(summary.loc[["Mean"]].round(6))
+    if args.full:
+        mean_row = summary.loc["Mean"]
+        for col in R_MEMBERSHIP_SUMMARY_MEAN:
+            if col in mean_row.index:
+                ref = R_MEMBERSHIP_SUMMARY_MEAN[col]
+                got = float(mean_row[col])
+                print(f"  {col}: got={got:.6f}, R={ref:.6f}, diff={abs(got - ref):.6f}")
+
+    print("\n[>] Property-based clustering (state + duration, max_clusters=5)...")
+    t3 = time.perf_counter()
+    pclust = property_based_clustering(
+        seqdata,
+        diss=diss,
+        properties=["state", "duration"],
+        max_clusters=5,
+        verbose=True,
+    )
+    print(f"    done in {time.perf_counter() - t3:.1f}s")
+    print_property_tree(pclust)
+
+    adj = pclust["info"].get("adjustment") or {}
+    r2 = adj.get("R2")
+    if r2 is not None:
+        print(f"\nGlobal R2: {r2:.5f}", end="")
+        if args.full:
+            print(f"  (R reference: {R_PROPERTY_GLOBAL_R2:.5f})")
+        else:
+            print()
+
+    print("\n[>] Cluster quality (as.clustrange)...")
+    pclustqual = property_clustering_quality(pclust, diss=diss, n_clusters=5)
+    print(pclustqual.stats.round(2))
+    if args.full:
+        _compare_table("cluster quality", pclustqual.stats, R_CLUSTER_QUALITY, ["PBC", "HG", "R2"])
+
+    print(f"\n[>] Total elapsed: {time.perf_counter() - t0:.1f}s")
+    if not args.full:
+        print(
+            "[i] Subsample run — use --full for R parity checks on membership / R2 / quality."
+        )
+
+
+if __name__ == "__main__":
+    main()

developer/important-literature/2018-studer-property-based-fuzzy-clustering.md

@@ -72,7 +72,7 @@
     "cluster_labels_from_kmedoids_result": ("sequenzo.clustering", "cluster_labels_from_kmedoids_result"),
     "hard_classification_variables": ("sequenzo.clustering", "hard_classification_variables"),
     "fanny_membership": ("sequenzo.clustering", "fanny_membership"),
-    "representative_indices_from_membership": ("sequenzo.clustering", "representative_indices_from_membership"),
+    "highest_membership_indices_from_membership": ("sequenzo.clustering", "highest_membership_indices_from_membership"),
     "soft_classification_variables": ("sequenzo.clustering", "soft_classification_variables"),
     "pseudoclass_regression": ("sequenzo.clustering", "pseudoclass_regression"),
     # multidomain
@@ -346,7 +346,7 @@ def __getattr__(name: str) -> Any:
     "cluster_labels_from_kmedoids_result",
     "hard_classification_variables",
     "fanny_membership",
-    "representative_indices_from_membership",
+    "highest_membership_indices_from_membership",
     "soft_classification_variables",
     "pseudoclass_regression",
     "create_idcd_sequence_from_csvs",

sequenzo/__init__.py

@@ -22,7 +22,8 @@
 from sequenzo.big_data.clara.utils.davies_bouldin import *
 from sequenzo.big_data.clara.utils.get_weighted_diss import get_weighted_diss
 from scipy.cluster.hierarchy import cut_tree
-from sequenzo.clustering.fuzzy.wfcmdd_fuzzy_clustering import wfcmdd
+
+from sequenzo.clustering.fuzzy_clustering import wfcmdd
 from sequenzo.clustering.k_medoids import KMedoids
 from sequenzo.clustering.sequences_to_variables.helske_regression_variables import (
     fanny_membership,

sequenzo/big_data/clara/clara.py

@@ -13,9 +13,15 @@ sequenzo/clustering/
   sequences_to_variables/
     helske_regression_variables.py # Helske (2024) sequence -> regression covariates
     helpers.py                     # shared small helpers
-  fuzzy/
+  fuzzy_clustering/
+    src/
+      fanny.cpp                    # R cluster::fanny (C++ port)
+      fanny.h
     wfcmdd_fuzzy_clustering.py     # distance-based fuzzy C-medoids
     fuzzy_sequence_plots.py        # membership-weighted seq index plots
+    fuzzy_helpers.py               # FANNY / wfcmdd entry points
+    fuzzy_regression.py            # Dirichlet / beta regression on membership
+  property_based_clustering/       # Studer (2018) seqpropclust workflow
   validation/
     partition_quality.py           # fixed partitions + CQI table
     dissmfacw_factors.py           # multi-factor discrepancy association core
@@ -41,9 +47,9 @@ sequenzo/clustering/
 | WeightedCluster `clustassoc` | `validation.cluster_covariate_association` | `cluster_association` | Clustering vs covariate |
 | WeightedCluster `rarcat` | `validation.rarcat_typology_regression` | `rarcat` | Robust typology AME |
 | `cluster::fanny` (membership) | `sequences_to_variables.helske_regression_variables` | `fanny_membership` | Distance-based soft membership |
-| WeightedCluster `wfcmdd` | `fuzzy.wfcmdd_fuzzy_clustering` | `wfcmdd` | FCMdd / NCdd / PCMdd |
-| WeightedCluster `fuzzyseqplot` | `fuzzy.fuzzy_sequence_plots` | `fuzzy_sequence_plot` | Membership-weighted index plot |
-| WeightedCluster `crispness` | `fuzzy.wfcmdd_fuzzy_clustering` | `crispness` | Partition sharpness |
+| WeightedCluster `wfcmdd` | `fuzzy_clustering.wfcmdd_fuzzy_clustering` | `wfcmdd` | FCMdd / NCdd / PCMdd |
+| WeightedCluster `fuzzyseqplot` | `fuzzy_clustering.fuzzy_sequence_plots` | `fuzzy_sequence_plot` | Membership-weighted index plot |
+| WeightedCluster `crispness` | `fuzzy_clustering.wfcmdd_fuzzy_clustering` | `crispness` | Partition sharpness |
 | Helske representativeness | `sequences_to_variables.helske_regression_variables` | `representativeness_matrix` | Not a TraMineR export |
 | Helske hard / soft / pseudoclass | `sequences_to_variables.helske_regression_variables` | `hard_classification_variables`, `soft_classification_variables`, `pseudoclass_regression` | Regression-ready typology covariates |
 

sequenzo/clustering/README.md

@@ -18,18 +18,40 @@
     FannyResult,
     fanny_membership,
     medoid_membership_approximation,
-    representative_indices_from_membership,
+    highest_membership_indices_from_membership,
     soft_classification_variables,
     pseudoclass_regression,
     max_distance,
     cluster_labels_to_dummies,
 )
-from .fuzzy import (
+from .fuzzy_clustering import (
     wfcmdd,
     crispness,
     WfcmddResult,
     fuzzy_sequence_plot,
     fuzzy_sequence_plot_single,
+    get_fuzzy_clusters,
+    FuzzyClusterResult,
+    membership_summary,
+    most_typical_members,
+    prepare_dirichlet_data,
+    DirichletRegData,
+    DirichletRegResult,
+    dirichlet_regression,
+    beta_regression,
+)
+from .property_based_clustering import (
+    extract_sequence_properties,
+    property_based_clustering,
+    seqpropclust,
+    cluster_split_schedule,
+    cut_tree,
+    prune_property_tree,
+    tree_labels,
+    property_clustering_quality,
+    print_property_tree,
+    plot_property_tree,
+    SUPPORTED_PROPERTIES,
 )
 from .validation import (
     cluster_range_from_partitions,
@@ -77,14 +99,34 @@ def _import_c_code():
     "FannyResult",
     "fanny_membership",
     "medoid_membership_approximation",
-    "representative_indices_from_membership",
+    "highest_membership_indices_from_membership",
     "soft_classification_variables",
     "pseudoclass_regression",
     "wfcmdd",
     "crispness",
     "WfcmddResult",
     "fuzzy_sequence_plot",
     "fuzzy_sequence_plot_single",
+    "get_fuzzy_clusters",
+    "FuzzyClusterResult",
+    "membership_summary",
+    "most_typical_members",
+    "prepare_dirichlet_data",
+    "DirichletRegData",
+    "DirichletRegResult",
+    "dirichlet_regression",
+    "beta_regression",
+    "extract_sequence_properties",
+    "property_based_clustering",
+    "seqpropclust",
+    "cluster_split_schedule",
+    "cut_tree",
+    "prune_property_tree",
+    "tree_labels",
+    "property_clustering_quality",
+    "print_property_tree",
+    "plot_property_tree",
+    "SUPPORTED_PROPERTIES",
     "cluster_range_from_partitions",
     "compute_partition_quality",
     "ClusterRangeResult",