Harmonic Sequence

biotuner/harmonic_sequence.py

Temporal harmonic-structure learning from successive biotuner outputs.

Module type: Objects

Six complementary approaches build on the outputs of compute_biotuner to learn or model the way harmonic content evolves over time:

1. HarmonicMarkov – Discrete-state Markov chain / HMM

2. WassersteinTrajectory – Optimal-transport trajectory in interval space

3. HarmonicDMD – Linear dynamical modes of metric vectors

4. HarmonicLatentSpace – Smooth manifold of harmonic states (PCA) 6. HarmonicTopology – Topological shape of harmonic trajectories (TDA) 7. HarmonicGrammar – N-gram / symbolic sequence model

All approaches are exposed through the orchestrating class HarmonicSequenceAnalyzer, which extracts representations from a list of compute_biotuner objects and feeds them to each model.

Quick start

>>> from biotuner.harmonic_sequence import HarmonicSequenceAnalyzer
>>> analyzer = HarmonicSequenceAnalyzer.from_biotuner_list(
...     bt_list, tuning="peaks_ratios"
... )
>>> analyzer.fit_all()
>>> print(analyzer.summary())

Module-level encoding functions (encode_histograms, encode_scalar_metrics, encode_ji_matrix) can also be imported and used standalone.

Optional dependencies

• scikit-learnKMeans, PCA, MDS, StandardScaler – required for

  HarmonicMarkov, HarmonicLatentSpace, WassersteinTrajectory.embed(), and HarmonicDMD.fit(use_histograms=True).

• hmmlearn : Gaussian HMM backend for HarmonicMarkov (use_hmm=True).

• ripserVietoris-Rips persistent homology for HarmonicTopology;

  gracefully falls back to H0-only via scipy linkage.

TUNING_ATTRS: Dict[str, str] = {'HE_scale': 'HE_scale', 'diss_scale': 'diss_scale', 'euler_fokker': 'euler_fokker', 'extended_peaks_ratios': 'extended_peaks_ratios', 'harm_fit_tuning': 'harm_fit_tuning_scale', 'harm_tuning': 'harm_tuning_scale', 'peaks_ratios': 'peaks_ratios', 'peaks_ratios_cons': 'peaks_ratios_cons'}

Maps user-facing tuning names to attribute names on compute_biotuner.

extract_tuning(bt_obj: Any, tuning: str = 'peaks_ratios') → List[float]

Extract a named tuning from a fitted compute_biotuner object.

Parameters:

• bt_obj (compute_biotuner) – A biotuner object after peaks_extraction() (and optionally compute_diss_curve(), compute_harmonic_entropy(), etc.).

• tuning (str) – One of the keys in TUNING_ATTRS.

Returns:

list of float – Positive finite ratios, or an empty list if the attribute is absent.

encode_histograms(ratios_list: List[List[float]], n_bins: int = 240, min_cents: float = 0.0, max_cents: float = 1200.0) → ndarray

Convert a sequence of ratio sets to a cents-histogram matrix.

Parameters:

• ratios_list (list of list of float, length T)

• n_bins (int, default=240) – Number of bins (5 cents/bin at default range).

• min_cents, max_cents (float) – Histogram range in cents.

Returns:

X (ndarray, shape (T, n_bins)) – Rows normalised to sum to 1 (zero-rows for empty ratio sets).

encode_scalar_metrics(bt_list: List[Any]) → ndarray

Build a scalar feature matrix from a list of compute_biotuner objects.

Extracts and concatenates:

• peaks_metrics: cons, tenney, harmsim, harm_fit, subharm_tension

• scale_metrics: diss_euler, dissonance, diss_harm_sim, diss_n_steps, HE, HE_n_steps, HE_harm_sim

• Mean peak frequency and peak count.

Parameters:

bt_list (list of compute_biotuner)

Returns:

X (ndarray, shape (T, D)) – Contains NaN wherever an attribute was absent or not computed.

clear_harmonicity_cache() → None

Empty the module-level dyad-matrix cache (forces recomputation).

encode_harmonicity_spectrum(bt_list: List[Any], *, fmin: float = 1.0, fmax: float = 30.0, precision_hz: float = 0.5, metric: str = 'harmsim', n_harms: int = 10, smoothness: int = 1, smoothness_harm: int = 1, normalize: bool = True, power_law_remove: bool = False, fs: float | None = None, cache: bool = True, progress: bool = False) → Tuple[ndarray, ndarray]

Compute the 1-D harmonicity spectrum per window (fast path).

For each compute_biotuner object, computes the per-frequency power-weighted dyad-similarity values. This is much faster than the full compute_global_harmonicity() because it:

• reuses the freq-only dyad-similarity matrix across all windows (it depends only on the frequency axis),

• skips phase, PLV, resonance, peak-finding, entropy, and DataFrame overhead (none of which we need here),

• caches results on each bt object under bt._harm_cache.

Parameters:

• bt_list (list of compute_biotuner) – Each bt must have bt.data (the raw signal) and bt.sf.

• fmin, fmax (float analysis range in Hz.)

• precision_hz (float frequency resolution; F = (fmax-fmin)/precision_hz.)

• metric (str ‘harmsim’ or ‘subharm_tension’.)

• n_harms (int number of harmonics in dyad similarity.)

• smoothness, smoothness_harm (int PSD / harmonicity-curve smoothing.)

• normalize (bool normalise by total power.)

• power_law_remove (bool remove 1/f trend before computation.)

• fs (float, optional override sampling frequency.)

• cache (bool, default=True reuse bt._harm_cache when present.)

• progress (bool print one line per window.)

Returns:

• H_spec (ndarray, shape (T, F))

• freqs (ndarray, shape (F,))

encode_harmonicity_matrices(bt_list: List[Any], *, fmin: float = 1.0, fmax: float = 30.0, precision_hz: float = 0.5, metric: str = 'harmsim', n_harms: int = 10, smoothness: int = 1, normalize: bool = True, power_law_remove: bool = False, fs: float | None = None, cache: bool = True, progress: bool = False) → Tuple[ndarray, ndarray]

Compute the 2-D harmonicity matrix per window (fast path).

Like encode_harmonicity_spectrum() but returns the F×F pairwise matrix per window — strictly more information. Shares the same cache and the same vectorised fast path.

Returns:

• M (ndarray, shape (T, F, F))

• freqs (ndarray, shape (F,))

encode_ji_matrix(ratios_list: List[List[float]], tolerance_cents: float = 30.0) → Tuple[ndarray, List[str]]

Binary-encode ratio sets as JI interval presence/absence vectors.

Each row is a binary vector over the interval_names dictionary: element j is 1 if any ratio in that timepoint’s set falls within tolerance_cents of interval j.

Parameters:

• ratios_list (list of list of float, length T)

• tolerance_cents (float, default=30) – Match radius in cents.

Returns:

• X (ndarray, shape (T, n_intervals) int8)

• interval_labels (list of str (column names))

histogram_to_ratios(h: ndarray, n_peaks: int | None = None, *, prominence: float = 0.0, min_cents: float = 0.0, max_cents: float = 1200.0, min_distance_cents: float = 25.0, include_unison: bool = True, include_octave: bool = True) → List[float]

Decode a cents histogram back to a list of frequency ratios.

Inverts encode_histograms() at the per-frame level by locating peaks in h and mapping their bin centres back to ratios via r = 2 ** (cents / 1200).

The same function works for observed histograms (one row of analyzer.histograms) and for synthesised histograms produced by Wasserstein barycenters, latent-space decode, DMD reconstruction, or Markov-state centroids.

Parameters:

• h (ndarray, shape (n_bins,)) – Normalised cents histogram. Sum should be 1; all-zero rows yield an empty (or unison/octave-only) ratio list.

• n_peaks (int or None, default=None) – If given, return at most this many peaks (top-amplitude). When None all local maxima above prominence are returned.

• prominence (float, default=0.0) – Minimum peak prominence (passed to scipy.signal.find_peaks). Useful to filter out noise in dense histograms (e.g. decoded from PCA latent space).

• min_cents, max_cents (float) – Histogram extent in cents (must match the encoder).

• min_distance_cents (float, default=25.0) – Minimum cents separation between detected peaks. Prevents adjacent bins from both registering as peaks.

• include_unison (bool, default=True) – Prepend 1.0 to the output (Scala convention; required for .scl).

• include_octave (bool, default=True) – Append 2.0 to the output (period of repetition).

Returns:

ratios (list of float) – Sorted ascending, finite, positive. [1.0] when h is all-zero and include_unison=True.

Examples

>>> # Round-trip: encode then decode a known tuning
>>> from biotuner.harmonic_sequence import (
...     encode_histograms, histogram_to_ratios,
... )
>>> h = encode_histograms([[1.5, 1.25, 1.333]])[0]
>>> histogram_to_ratios(h, include_unison=False, include_octave=False)
[1.2497..., 1.3325..., 1.4983...]

histogram_to_scl(h: ndarray, name: str = 'biotuner_histogram', *, n_peaks: int | None = None, prominence: float = 0.0, write: bool = False, **kwargs: Any) → str

Render a cents histogram as a Scala (.scl) tuning-file string.

Wraps histogram_to_ratios() and biotuner.biotuner_utils.create_SCL(). Pass write=True to also save <name>.scl to disk.

Parameters:

• h (ndarray, shape (n_bins,))

• name (str) – Used both as the in-file scale title and (when write=True) the filename stem.

• n_peaks, prominence (passed to histogram_to_ratios().)

• write (bool, default=False) – When True, write <name>.scl next to the working directory.

• **kwargs – Forwarded to histogram_to_ratios() (e.g. min_distance_cents).

Returns:

scl_str (str) – Scala-formatted text.

histograms_to_midi(H: ndarray, filename: str = 'harmonic_sequence', *, base_freq: float = 220.0, duration_beats: float = 1.0, velocity: int = 80, n_peaks: int | None = 5, prominence: float = 0.0, microtonal: bool = True, subdivision: int = 1, skip_empty: bool = True, **kwargs: Any) → Any

Render a sequence of histograms as a microtonal MIDI file.

Each row of H becomes one chord: peaks of the histogram are decoded to ratios, multiplied by base_freq to get absolute frequencies, then sent to biotuner.biotuner_utils.create_midi() (one channel per voice with pitch bends for microtonal precision).

The function works equally well on:

• analyzer.histograms — the recorded biosignal trajectory

• analyzer.wasserstein.interpolate_pair(t1, t2, n) — a glissando

• analyzer.latent.decode(Z_path) — a synthesised path through harmonic space

• analyzer.markov._km.cluster_centers_ — prototype tunings

Parameters:

• H (ndarray, shape (T, n_bins) or (n_bins,)) – Sequence of histograms. A single 1-D histogram is treated as one chord. Rows do not need to be normalised.

• filename (str, default=’harmonic_sequence’) – Output filename stem; .mid is appended by create_midi.

• base_freq (float, default=220.0) – Frequency that the unison ratio (1.0) maps to (Hz). 220 Hz = A3.

• duration_beats (float, default=1.0) – Note duration per chord. Pass a list/array of length T for variable durations (e.g. modulated by Wasserstein flux).

• velocity (int, default=80) – MIDI velocity (1-127) for every note.

• n_peaks (int or None, default=5) – Maximum peaks decoded per histogram. None keeps all peaks.

• prominence (float, default=0.0) – Minimum peak prominence (helps for smooth/decoded histograms).

• microtonal (bool, default=True) – Emit pitch-bend messages so non-12-TET ratios play in tune.

• subdivision (int, default=1) – Beat subdivisions (passed to create_midi).

• skip_empty (bool, default=True) – Drop windows whose histogram has no detected peaks. When False an empty chord is inserted (silent rest of duration_beats).

• **kwargs – Forwarded to histogram_to_ratios().

Returns:

mid (mido.MidiFile) – The saved MIDI file object.

Notes

With microtonal=True, each voice occupies its own MIDI channel, capping the chord size at 16 notes. n_peaks=5 is a safe default.

find_optimal_n_states(X: ndarray, k_range: Tuple[int, int] = (2, 10), random_state: int = 42) → Tuple[int, Dict[int, float]]

Select the number of Markov states via silhouette score.

For each K in k_range, K-means is fitted and the mean silhouette coefficient is computed. The K with the highest score is returned.

Parameters:

• X (ndarray, shape (T, n_features)) – Histogram (or other feature) matrix.

• k_range ((int, int), default=(2, 10)) – Inclusive range of candidate K values.

• random_state (int, default=42)

Returns:

• best_k (int)

• scores (dict {k: silhouette_score})

Notes

Requires scikit-learn. Silhouette ranges from -1 (poor) to +1 (ideal); values above ~0.5 indicate well-separated clusters.

class HarmonicMarkov(n_states: int = 5, order: int = 1, auto_k_range: Tuple[int, int] = (2, 10), use_hmm: bool = False, random_state: int = 42)

Bases: object

Discrete-state Markov chain (and optional HMM) over harmonic states.

Each timepoint is assigned to one of n_states clusters defined by KMeans on the cents-histogram representation. A Markov chain of arbitrary order is estimated from the state sequence.

Parameters:

• n_states (int or 'auto', default=5) – Number of discrete harmonic states. When 'auto', the optimal K is selected automatically via silhouette score over the range given by auto_k_range.

• order (int, default=1) – Markov chain memory. order=1 is a standard first-order chain (next state depends only on the current state). order=2 conditions on the last two states, etc. Higher orders can capture longer harmonic patterns but require more data (at least n_states**order transitions are needed for reliable estimates).

• auto_k_range ((int, int), default=(2, 10)) – Range of K to search when n_states='auto'.

• use_hmm (bool, default=False) – Fit a Gaussian HMM (requires hmmlearn).

• random_state (int, default=42)

fit(X: ndarray) → HarmonicMarkov

Fit the Markov chain on a histogram matrix.

Parameters:

X (ndarray, shape (T, n_bins))

Returns:

self

property transition_matrix_: ndarray

Row-normalised (n_states × n_states) first-order transition matrix.

property high_order_transition_: Dict[tuple, ndarray] | None

Higher-order transition table.

None when order=1. Otherwise a dict mapping each observed history tuple of length order to a probability distribution (ndarray of shape (n_states,)) over the next state.

property state_labels_: ndarray

Integer state label for each fitted timepoint.

property silhouette_scores_: Dict[int, float] | None

Silhouette scores per K when n_states='auto', else None.

property steady_state_: ndarray

Stationary distribution π such that πT = π (first-order).

property transition_entropy_: float

Mean Shannon entropy (bits) of the transition distributions.

For order=1 this is the mean row entropy of the K×K matrix. For higher orders it is the mean entropy over all observed history tuples, weighted equally.

predict_next_proba(history) → ndarray

Transition probability distribution from the given history.

Parameters:

history (int or tuple of int) – For order=1, a single state index. For order>1, a tuple of the last order state indices (oldest first).

Returns:

proba (ndarray, shape (n_states,)) – Falls back to the first-order row when the history tuple was not observed during training.

decode_viterbi(X: ndarray) → ndarray

Viterbi decoding via HMM; falls back to KMeans if HMM unavailable.

class WassersteinTrajectory(n_bins: int = 240)

Bases: object

Optimal-transport trajectory in interval space.

Treats each timepoint’s cents histogram as a probability distribution over [0, 1200] and computes the 1-D Wasserstein (Earth Mover’s) distance between frames. The pairwise distance matrix supports manifold visualisation; consecutive distances measure harmonic flux.

Parameters:

n_bins (int, default=240) – Number of histogram bins (must match the encoder).

fit(X: ndarray) → WassersteinTrajectory

Compute pairwise Wasserstein distances and per-frame flux.

Parameters:

X (ndarray, shape (T, n_bins)) – Normalised cents histograms (rows should sum to 1).

Returns:

self

property distance_matrix_: ndarray

Pairwise Wasserstein distance matrix (T × T).

property flux_: ndarray

Per-frame harmonic velocity (consecutive W1 distances), length T-1.

barycenter(h1: ndarray, h2: ndarray, alpha: float) → ndarray

1-D Wasserstein barycenter between h1 and h2 at weight alpha.

Interpolates in quantile space: Q_α = (1-α)·Q₁ + α·Q₂, then maps back to a histogram via sampling.

Parameters:

• h1, h2 (ndarray, shape (n_bins,) normalised histograms)

• alpha (float in [0, 1])

Returns:

bary (ndarray, shape (n_bins,) normalised)

interpolate_pair(t1: int, t2: int, n_steps: int = 10) → List[ndarray]

Return n_steps histograms interpolating between timepoints t1 and t2.

Parameters:

• t1, t2 (int indices into the fitted sequence)

• n_steps (int)

Returns:

list of ndarray, each shape (n_bins,)

embed(n_components: int = 2, method: str = 'mds') → ndarray

Low-dimensional embedding of the Wasserstein distance matrix.

Parameters:

• n_components (int, default=2)

• method (str, default=’mds’) – 'mds' (requires sklearn) or 'umap' (requires umap-learn).

Returns:

Z (ndarray, shape (T, n_components))

class HarmonicDMD(rank: int | None = None, center: bool = True)

Bases: object

Linear dynamical modes of harmonic metric vectors.

Applies the exact DMD algorithm to a (T × D) feature matrix. Eigenvalues inside the unit circle correspond to decaying transients; on the unit circle to oscillatory harmonic cycles; outside to growing instabilities.

Parameters:

• rank (int or None, default=None) – SVD truncation rank. None keeps all singular values above a numerical threshold (relative to the largest singular value).

• center (bool, default=True) – Subtract the temporal mean before fitting.

fit(X: ndarray) → HarmonicDMD

Fit the DMD model.

Parameters:

X (ndarray, shape (T, D)) – Feature matrix. NaN values are imputed with column means.

Returns:

self

property eigenvalues_: ndarray

DMD eigenvalues (complex array, length r).

property modes_: ndarray

DMD spatial modes (D × r, complex).

property growth_rates_: ndarray

positive → growing, negative → decaying.

Type:

Real part of log(λ)

property frequencies_: ndarray

Imag part of log(λ), proportional to oscillation frequency.

oscillatory_modes(threshold: float = 0.05) → Tuple[ndarray, ndarray]

Return eigenvalues and indices of modes close to the unit circle.

Parameters:

threshold (float, default=0.05) – Maximum allowed distance from |λ| = 1.

Returns:

• eigenvalues (complex ndarray)

• idx (int ndarray)

reconstruct(n_steps: int) → ndarray

Forward-propagate from the last fitted state for n_steps.

Uses the DMD mode expansion: x(k) = Φ · diag(λ)^k · b₀ where b₀ = Φ⁺ x₀.

Returns:

X_pred (ndarray, shape (n_steps, D) (with mean added back))

class HarmonicLatentSpace(latent_dim: int = 3, random_state: int = 42)

Bases: object

Smooth low-dimensional manifold of harmonic states via PCA.

Fits a linear PCA to the cents-histogram matrix, providing:

• encode(X) – project histograms into latent coordinates

• decode(Z) – reconstruct histograms from latent coordinates

• interpolate(z1, z2) – smooth path between two harmonic states

Parameters:

• latent_dim (int, default=3)

• random_state (int, default=42)

fit(X: ndarray) → HarmonicLatentSpace

Fit PCA to the histogram matrix.

Parameters:

X (ndarray, shape (T, n_bins))

Returns:

self

encode(X: ndarray) → ndarray

Project histograms into latent space.

Parameters:

X (ndarray, shape (T, n_bins))

Returns:

Z (ndarray, shape (T, latent_dim))

decode(Z: ndarray) → ndarray

Reconstruct histograms from latent coordinates.

Parameters:

Z (ndarray, shape (T, latent_dim))

Returns:

X_rec (ndarray, shape (T, n_bins) (clipped to non-negative))

trajectory() → ndarray

Latent coordinates of all fitted timepoints, shape (T, latent_dim).

interpolate(z1: ndarray, z2: ndarray, n_steps: int = 10) → ndarray

Linear interpolation between two latent points.

Parameters:

• z1, z2 (ndarray, shape (latent_dim,))

• n_steps (int)

Returns:

Z_path (ndarray, shape (n_steps, latent_dim))

property explained_variance_ratio_: ndarray

Fraction of variance explained by each latent dimension.

property reconstruction_error_: float

Mean squared reconstruction error on the fitted data.

class HarmonicTopology(embedding_dim: int = 3, delay: int = 1)

Bases: object

Topological shape of harmonic scalar trajectories.

Applies Takens delay embedding to convert a scalar harmonicity time series (e.g., harmsim, subharm_tension) into a point cloud, then computes persistent homology:

• H0 (connected components): harmonic attractor basins

• H1 (loops): cyclic / periodic harmonic progressions

H0 is always available via scipy single-linkage clustering. H1 requires ripser (pip install ripser).

Parameters:

• embedding_dim (int, default=3) – Takens embedding dimension d.

• delay (int, default=1) – Takens delay τ (in time-steps).

fit(scalar_series: ndarray) → HarmonicTopology

Embed the scalar series and compute persistent homology.

Parameters:

scalar_series (array-like, shape (T,)) – A 1-D harmonicity time series. NaN values are linearly interpolated.

Returns:

self

property takens_embedding_: ndarray

Point cloud from Takens delay embedding, shape (N, embedding_dim).

property persistence_diagram_: List[ndarray]

List of persistence diagrams [H0, H1] (H1 only if ripser installed).

Each element is an ndarray of (birth, death) pairs.

property betti_numbers_: ndarray

Betti numbers β0, β1 evaluated at the median finite filtration value.

Returns:

betti (ndarray of int, shape (2,))

session_fingerprint() → ndarray

Fixed-length descriptor summarising the persistence diagrams.

Returns a 6-element vector: [mean_H0_pers, max_H0_pers, n_H0_bars,

mean_H1_pers, max_H1_pers, n_H1_bars]

class HarmonicGrammar(tolerance_cents: float = 30.0, n_gram: int = 2)

Bases: object

N-gram symbolic sequence model over JI interval-labelled chords.

Each timepoint’s ratio set is mapped to a frozenset of JI interval names (using the biotuner interval_names dictionary), creating a symbolic chord sequence. An N-gram model captures the conditional probability of each chord token given the preceding N-1 tokens.

Parameters:

• tolerance_cents (float, default=30.0) – Match radius for JI interval labelling.

• n_gram (int, default=2) – N-gram order (2 = bigram, 3 = trigram, …).

fit(ratios_list: List[List[float]]) → HarmonicGrammar

Label each ratio set and build the N-gram language model.

Parameters:

ratios_list (list of list of float, length T)

Returns:

self

property chord_sequence_: List[frozenset]

Symbolic chord sequence (list of frozensets of interval names).

property vocabulary_: List[frozenset]

Unique chord tokens in order of first occurrence.

transition_proba(context: Tuple) → Dict[frozenset, float]

Conditional probability P(next_chord | context).

Parameters:

context (tuple of frozenset(s)) – The preceding (n_gram - 1) chords.

Returns:

dict mapping frozenset → probability

property transition_entropy_: float

Mean Shannon entropy of (n-1)-gram conditional distributions (bits).

top_ngrams(top_k: int = 10) → List[Tuple[Tuple, int]]

Return the top_k most frequent N-grams with their counts.

top_motifs(min_length: int = 2, max_length: int = 4, top_k: int = 10) → List[Tuple[Tuple, int]]

Find the most frequent chord motifs of lengths in [min_length, max_length].

Returns:

list of (motif_tuple, count) sorted by frequency descending.

static levenshtein(seq1: List[frozenset], seq2: List[frozenset]) → int

Levenshtein edit distance between two symbolic chord sequences.

Substitution cost is 1 for non-identical chords (frozenset equality).

Parameters:

seq1, seq2 (list of frozenset)

Returns:

int

class HarmonicSequenceAnalyzer(n_hist_bins: int = 240, tolerance_cents: float = 30.0, representation: str = 'cents_histogram', representation_kwargs: Dict[str, Any] | None = None)

Bases: object

Orchestrating class for temporal harmonic-structure learning.

Extracts representations from a list of compute_biotuner objects (or raw ratio lists) and exposes each modelling approach via a consistent fit_* interface with lazy encoding caches.

Parameters:

• n_hist_bins (int, default=240) – Histogram resolution for cents encoding (5 cents/bin over one octave).

• tolerance_cents (float, default=30.0) – JI matching tolerance for grammar and JI matrix encoding.

Examples

>>> analyzer = HarmonicSequenceAnalyzer.from_biotuner_list(
...     bt_list, tuning="peaks_ratios"
... )
>>> analyzer.fit_all()
>>> print(analyzer.summary())

VALID_REPRESENTATIONS = ('cents_histogram', 'harmonicity_spectrum', 'harmonicity_matrix')

classmethod from_biotuner_list(bt_list: List[Any], tuning: str = 'peaks_ratios', n_hist_bins: int = 240, tolerance_cents: float = 30.0, representation: str = 'cents_histogram', representation_kwargs: Dict[str, Any] | None = None) → HarmonicSequenceAnalyzer

Build an analyzer from a list of fitted compute_biotuner objects.

Parameters:

• bt_list (list of compute_biotuner) – Each element should have had at least peaks_extraction() called. Additional attributes (peaks_metrics, scale_metrics, diss_scale, etc.) are used when present.

• tuning (str, default=’peaks_ratios’) – Which tuning attribute to use as the primary ratio source. Supported: 'peaks_ratios', 'peaks_ratios_cons', 'diss_scale', 'HE_scale', 'euler_fokker', 'harm_tuning', 'harm_fit_tuning', 'extended_peaks_ratios'.

• n_hist_bins (int, default=240)

• tolerance_cents (float, default=30.0)

• representation (str, default=’cents_histogram’) – Feature representation consumed by Markov / Wasserstein / DMD / Latent / Topology. One of:

  • 'cents_histogram' : 240-bin cents distribution (default).

  • 'harmonicity_spectrum' : 1-D power-weighted harmonicity per frequency. Requires raw signal in bt.data. Preserves frequency identity (alpha vs gamma).

  • 'harmonicity_matrix' : 2-D F×F pairwise harmonicity matrix flattened to F² per window. Strictly more information than the spectrum (no collapse along the second frequency axis).

  Grammar is always cents-based (it labels ratios via JI dictionary). The MIDI/SCL rendering bridge always uses cents histograms, regardless of this parameter.

• representation_kwargs (dict, optional) – Forwarded to encode_harmonicity_spectrum() / encode_harmonicity_matrices(). Common keys: fmin, fmax, precision_hz, metric, n_harms.

Returns:

HarmonicSequenceAnalyzer

classmethod from_ratios_list(ratios_list: List[List[float]], n_hist_bins: int = 240, tolerance_cents: float = 30.0) → HarmonicSequenceAnalyzer

Build an analyzer directly from a list of ratio sets (no bt objects).

Parameters:

• ratios_list (list of list of float, length T)

• n_hist_bins (int, default=240)

• tolerance_cents (float, default=30.0)

Returns:

HarmonicSequenceAnalyzer

property histograms: ndarray

Cents-histogram matrix (T × n_hist_bins), computed on first access.

property scalar_metrics: ndarray

Scalar metric matrix (T × D) from bt objects, or empty array.

property features: ndarray

Active feature matrix (T, D) consumed by the fit methods.

Determined by self.representation:

• 'cents_histogram' → histograms (T × n_hist_bins).

• 'harmonicity_spectrum' → (T × F) harmonicity per frequency.

• 'harmonicity_matrix' → (T × F²) flattened pairwise matrices.

Cached on first access. For the unflattened (T × F × F) tensor of the matrix representation, see harmonicity_matrices_.

property harmonicity_matrices_: ndarray | None

The (T, F, F) tensor when representation='harmonicity_matrix'.

None for other representations.

property features_freqs_: ndarray | None

Frequency axis (Hz) for the harmonicity representations, else None.

fit_markov(n_states=5, order: int = 1, auto_k_range: Tuple[int, int] = (2, 10), use_hmm: bool = False, random_state: int = 42, **kwargs) → HarmonicMarkov

Fit a Markov chain (and optionally HMM) on the active feature matrix.

Consumes features — i.e. cents histograms, the harmonicity spectrum, or the flattened harmonicity matrix, depending on self.representation.

Parameters:

• n_states (int or 'auto', default=5) – Number of states. Pass 'auto' to select K automatically via silhouette score (see find_optimal_n_states()).

• order (int, default=1) – Markov chain order (memory length). order=1 is the standard first-order chain; higher values condition on the last order states.

• auto_k_range ((int, int), default=(2, 10)) – Search range when n_states='auto'.

• use_hmm (bool, default=False)

• random_state (int, default=42)

Returns:

HarmonicMarkov (also stored as self.markov)

fit_wasserstein() → WassersteinTrajectory

Compute pairwise Wasserstein distances on the active features.

For the cents-histogram representation, distances are over the cents axis (interval space). For the harmonicity representations, each row is first L1-normalised so it is treated as a probability density over frequency (or pairwise frequency for the matrix mode).

Returns:

WassersteinTrajectory (also stored as self.wasserstein)

fit_dmd(rank: int | None = None, center: bool = True, use_histograms: bool = False, **kwargs) → HarmonicDMD

Fit DMD on the scalar metric (or histogram PCA) feature matrix.

Parameters:

• rank (int or None) – SVD truncation rank.

• center (bool, default=True) – Subtract temporal mean before fitting.

• use_histograms (bool, default=False) – If True, apply PCA (10 components) to the histograms first and run DMD on the reduced representation. Requires sklearn.

Returns:

HarmonicDMD (also stored as self.dmd)

fit_latent(latent_dim: int = 3, random_state: int = 42, **kwargs) → HarmonicLatentSpace

Fit a PCA-based latent harmonic space on the histogram encoding.

Parameters:

• latent_dim (int, default=3)

• random_state (int, default=42)

Returns:

HarmonicLatentSpace (also stored as self.latent)

fit_topology(scalar_key: str = 'harmsim', embedding_dim: int = 3, delay: int = 1, **kwargs) → HarmonicTopology

Fit TDA on a scalar harmonicity time series.

Parameters:

• scalar_key (str, default=’harmsim’) – Which scalar to use for the Takens embedding. Options:

  • Any key in peaks_metrics: 'harmsim', 'cons', 'tenney', 'harm_fit', 'subharm_tension'

  • Any key in scale_metrics: 'dissonance', 'HE', 'diss_harm_sim', etc.

  • 'latent_N' (e.g. 'latent_0') uses the N-th dimension of a previously fitted HarmonicLatentSpace model.

• embedding_dim (int, default=3)

• delay (int, default=1)

Returns:

HarmonicTopology (also stored as self.topology)

fit_grammar(n_gram: int = 2, **kwargs) → HarmonicGrammar

Fit an N-gram interval grammar on the ratio sequence.

Parameters:

n_gram (int, default=2 (bigram))

Returns:

HarmonicGrammar (also stored as self.grammar)

fit_all(n_states=5, markov_order: int = 1, auto_k_range: Tuple[int, int] = (2, 10), latent_dim: int = 3, n_gram: int = 2, topology_scalar: str = 'harmsim', **kwargs) → HarmonicSequenceAnalyzer

Run all six modelling approaches, silently skipping any that fail.

Parameters:

• n_states (int or 'auto', default=5) – Number of Markov states. 'auto' triggers automatic K selection via silhouette score over auto_k_range.

• markov_order (int, default=1) – Markov chain memory depth.

• auto_k_range ((int, int), default=(2, 10)) – K search range when n_states='auto'.

• latent_dim (int, default=3) – Latent space dimensionality.

• n_gram (int, default=2) – Grammar N-gram order.

• topology_scalar (str, default=’harmsim’) – Scalar key for the TDA embedding.

Returns:

self

get_histograms(source: str = 'observed', *, t1: int | None = None, t2: int | None = None, n_steps: int = 16, n_predict: int = 16) → ndarray

Pull a histogram sequence from any source — observed or generated.

Single entry point for the rendering bridge. Returns a (T, n_bins) matrix that can be fed to histograms_to_midi() / histogram_to_scl().

Parameters:

• source (str) – One of:

  • 'observed' : the recorded histograms (one per window)

  • 'wasserstein_interp' : Wasserstein barycenter path between t1 and t2 (requires fit_wasserstein())

  • 'latent_interp' : decoded path between latent points at t1 and t2 (requires fit_latent())

  • 'markov_centroids' : K prototype histograms (KMeans cluster centres; requires fit_markov())

  • 'markov_sample' : sample a state sequence of length n_steps and return the corresponding centroids

  • 'dmd_predict' : forward-extrapolate n_predict steps (requires DMD fitted with use_histograms=True)

• t1, t2 (int, optional) – Window indices for interpolation modes (default: 0 and T-1).

• n_steps (int, default=16) – Length of the interpolation / Markov sample.

• n_predict (int, default=16) – Steps to extrapolate for 'dmd_predict'.

Returns:

H (ndarray, shape (T_out, n_hist_bins))

to_scl(index: int = 0, name: str | None = None, source: str = 'observed', *, n_peaks: int | None = None, prominence: float = 0.0, write: bool = False, **kwargs: Any) → str

Export one histogram of the sequence as a Scala (.scl) tuning.

Parameters:

• index (int, default=0) – Row of the histogram source to export.

• name (str, optional) – Scale name (also filename when write=True). Defaults to f"biotuner_{source}_{index}".

• source (str) – See get_histograms().

• n_peaks, prominence, write – Forwarded to histogram_to_scl().

Returns:

str Scala file contents.

to_midi(filename: str = 'harmonic_sequence', source: str = 'observed', *, base_freq: float = 220.0, duration_beats: float = 1.0, velocity: int = 80, n_peaks: int | None = 5, prominence: float = 0.0, flux_modulated_durations: bool = False, **kwargs: Any) → Any

Export a histogram sequence as a microtonal MIDI file.

Parameters:

• filename (str) – Output filename stem (.mid is appended).

• source (str) – See get_histograms().

• base_freq (float, default=220.0) – Frequency that ratio 1.0 maps to (Hz).

• duration_beats (float, default=1.0) – Beats per chord. Ignored when flux_modulated_durations=True.

• velocity (int, default=80)

• n_peaks (int or None, default=5) – Maximum simultaneous voices per chord (also the MIDI channel cap when microtonal=True).

• prominence (float, default=0.0)

• flux_modulated_durations (bool, default=False) – When True (and Wasserstein is fitted), each chord’s duration is inversely proportional to local flux: high-flux moments get shorter notes, stable moments get longer. Useful for sonifying the harmonic-velocity envelope alongside the tunings themselves.

• **kwargs – Forwarded to get_histograms() (e.g. t1, t2, n_steps).

Returns:

mido.MidiFile

summary() → str

Return a concise text summary of all fitted models.