Chickaboo/Advanced-MIDI-Renderer
1
1#! /usr/bin/python32 3r'''############################################################################4################################################################################5#6#7# Tegridy Cupy Python Module (TCUPY)8# Version 1.09#10# Project Los Angeles11#12# Tegridy Code 202613#14# https://github.com/asigalov61/tegridy-tools15#16#17################################################################################18#19# Copyright 2026 Project Los Angeles / Tegridy Code20#21# Licensed under the Apache License, Version 2.0 (the "License");22# you may not use this file except in compliance with the License.23# You may obtain a copy of the License at24#25# http://www.apache.org/licenses/LICENSE-2.026#27# Unless required by applicable law or agreed to in writing, software28# distributed under the License is distributed on an "AS IS" BASIS,29# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.30# See the License for the specific language governing permissions and31# limitations under the License.32#33################################################################################34################################################################################35#36# Critical dependencies37#38# !pip install cupy-cuda13x39# !pip install numpy==1.26.440#41################################################################################42'''43 44################################################################################45 46print('=' * 70)47print('Loading module...')48print('Please wait...')49print('=' * 70)50 51################################################################################52 53import sys54import os55import tqdm56 57################################################################################58 59try:60 import cupy as cp61 import numpy as np62 print('=' * 70)63 print('CuPy is found!')64 print('Will use CuPy and GPU for processing!')65 print('=' * 70)66 67except ImportError as e:68 print('Error: Could not import CuPy!')69 print(f'Details: {e}')70 print('=' * 70)71 print("Please make sure CuPy is installed.")72 print('pip install cupy-cuda13x')73 print('=' * 70)74 print('Will use NumPy for now...')75 import numpy as cp76 import numpy as np77 print('=' * 70)78 79################################################################################80 81from collections import defaultdict, deque82from typing import Optional, Tuple, Dict, Any, List83 84################################################################################85 86# Constants87MEMORY_LEN = 12 # Autoregressive context length88SEQUENCE_LENGTH = 32 # Each sequence has 24 triplets89 90# Baseline penalty values:91REPETITION_PENALTY = (1.0, 1.0, 1.0) # base repetition penalty per element92SPIKE_PENALTY_STRENGTH = (1.0, 1.0, 1.0) # base spike penalty strength per element93SPIKE_SIGMA = (1.0, 1.0, 1.0) # baseline sigma value per element (minimum allowed)94 95###################################################################################96 97def find_numpy_array(src_array, trg_array):98 99 """100 Finds 1D numpy array in 2D numpy array101 """102 103 match_mask = np.all(src_array == trg_array, axis=1)104 105 return np.where(match_mask)[0]106 107###################################################################################108 109def vertical_list_search(src_list, trg_list):110 111 """112 For each vertical window of consecutive rows of height len(trg_list) in src_list,113 this function checks whether for every offset j (0 <= j < len(trg_list)) the row114 at index (window_start + j) contains trg_list[j].115 116 It returns a list of windows (each a list of consecutive row indices) that meet this condition.117 """118 119 if not src_list or not trg_list:120 return []121 122 n = len(src_list)123 k = len(trg_list)124 125 num_windows = n - k + 1126 127 if num_windows <= 0:128 return []129 130 # Determine the maximum row length.131 max_len = max(len(row) for row in src_list)132 133 # Determine a fill value guaranteed to be less than any valid value.134 global_min = min(min(row) for row in src_list if row)135 fill_value = global_min - 1136 137 # Build a padded 2D array A (shape n x max_len) from src_list.138 A = np.full((n, max_len), fill_value, dtype=np.int64)139 for i, row in enumerate(src_list):140 L = len(row)141 A[i, :L] = row142 143 # For each unique target in trg_list, compute a Boolean vector of length n.144 # present[t][i] will be True if A[i, :] contains t, else False.145 unique_targets = set(trg_list)146 147 present_dict = {}148 149 for t in unique_targets:150 # Compute along axis=1 so that for each row we see if any element equals t.151 present_dict[t] = np.any(A == t, axis=1)152 153 # Build a Boolean array B of shape (k, num_windows) where for each offset j,154 # B[j, s] = present_dict[ trg_list[j] ][s + j] for each window starting index s.155 B = np.empty((k, num_windows), dtype=bool)156 157 for j in range(k):158 t = trg_list[j]159 # For a vertical window starting at s, row s+j should contain t.160 B[j, :] = present_dict[t][j: j + num_windows]161 162 # A window is valid if all k rows in that window contain the required target.163 valid_windows_mask = np.all(B, axis=0)164 valid_starts = np.nonzero(valid_windows_mask)[0]165 166 # Create output windows (each as a list of consecutive row indices).167 result = [list(range(s, s + k)) for s in valid_starts]168 169 return result170 171 172###################################################################################173 174def pack_sequences(train_data, pad_val=-1):175 """176 Packs a list of variable-length token sequences into a 2D CuPy array.177 178 This version computes lengths and builds the padded array and mask entirely on GPU.179 It converts each sequence into a CuPy array, concatenates them, and assigns tokens in one shot.180 181 Returns:182 batch: a CuPy array of shape (n, max_len)183 lengths: a CuPy array of shape (n,) containing each sequence's length.184 """185 n = len(train_data)186 # Compute lengths of each sequence and convert to a CuPy array.187 lengths = cp.array([len(seq) for seq in train_data], dtype=cp.int64)188 max_len_val = int(cp.max(lengths).get())189 # Allocate the padded 2D array filled with pad_val.190 batch = cp.full((n, max_len_val), pad_val, dtype=cp.int64)191 # Create a boolean mask: for each row, positions less than the sequence length are valid.192 mask = cp.arange(max_len_val).reshape(1, max_len_val) < lengths.reshape(n, 1)193 # Convert each sequence to a CuPy array and concatenate them.194 sequences = [cp.array(seq, dtype=cp.int64) for seq in train_data]195 flat = cp.concatenate(sequences)196 # Fill in the valid positions.197 batch[mask] = flat198 return batch, lengths199 200###################################################################################201 202def count_best_pair_gpu(batch, lengths, factor, pad_val=-1):203 """204 Given the entire GPU-resident packed data, compute the most frequent205 adjacent pair (encoded as: pair_val = first * factor + second) on GPU.206 """207 n, L = batch.shape208 cols = cp.arange(L - 1, dtype=cp.int64)209 cols_expanded = cp.broadcast_to(cols, (n, L - 1))210 valid_mask = cols_expanded < cp.reshape(lengths, (n, 1)) - 1211 212 first_tokens = batch[:, :L - 1]213 second_tokens = batch[:, 1:L]214 valid_first = first_tokens[valid_mask]215 valid_second = second_tokens[valid_mask]216 217 pairs = valid_first * factor + valid_second218 if pairs.size == 0:219 return None220 221 sorted_pairs = cp.sort(pairs)222 diff = cp.diff(sorted_pairs)223 boundaries = cp.nonzero(diff)[0] + 1224 group_starts = cp.concatenate([cp.array([0], dtype=cp.int64), boundaries])225 group_ends = cp.concatenate([boundaries, cp.array([sorted_pairs.size], dtype=cp.int64)])226 group_counts = group_ends - group_starts227 228 max_idx = int(cp.argmax(group_counts))229 best_pair_enc = int(sorted_pairs[group_starts[max_idx]])230 best_freq = int(group_counts[max_idx])231 first = best_pair_enc // factor232 second = best_pair_enc % factor233 return (first, second, best_freq)234 235###################################################################################236 237merge_kernel_code = r'''238extern "C" __global__239void merge_pair_kernel(const long* input, long* output, 240 const long* input_lengths, long* output_lengths,241 const long num_rows, const long num_cols,242 const long a, const long b, const long new_token,243 const long pad_val) {244 int row = blockIdx.x * blockDim.x + threadIdx.x;245 if (row >= num_rows) return;246 long in_length = input_lengths[row];247 long out_idx = 0;248 bool skip_next = false;249 for (long i = 0; i < in_length; i++) {250 if (skip_next) {251 skip_next = false;252 continue;253 }254 long token = input[row * num_cols + i];255 if (i < in_length - 1 && token == a && input[row * num_cols + i + 1] == b) {256 output[row * num_cols + out_idx] = new_token;257 out_idx++;258 skip_next = true;259 } else {260 output[row * num_cols + out_idx] = token;261 out_idx++;262 }263 }264 output_lengths[row] = out_idx;265 for (long j = out_idx; j < num_cols; j++) {266 output[row * num_cols + j] = pad_val;267 }268}269'''270 271try:272 merge_kernel = cp.RawKernel(merge_kernel_code, 'merge_pair_kernel')273 274except:275 pass276 277###################################################################################278 279def learn_bpe_codes_gpu(train_data, vocab_size=4096, max_merges=None, pad_val=-1):280 """281 Learn BPE merge rules completely on GPU.282 283 The training data is packed once (using the vectorized pack_sequences).284 On each merge iteration, the best adjacent pair is computed on GPU and then merged285 into a new token via a custom merge kernel (with double-buffering).286 287 Returns:288 codes: a list of merge rules as ((first, second), new_token)289 final_data: the merged training data (list of sequences)290 """291 # Pack the entire dataset onto GPU.292 batch, lengths = pack_sequences(train_data, pad_val)293 n, L = batch.shape294 295 # Initialize vocabulary and the next available token.296 initial_vocab = {token for seq in train_data for token in seq}297 next_token = max(initial_vocab) + 1298 codes = []299 merge_count = 0300 pbar = tqdm.tqdm(total=max_merges if max_merges is not None else None,301 desc="Learning BPE Codes (GPU)", leave=True)302 303 # Preallocate buffers for double-buffering.304 work_batch = cp.empty_like(batch)305 work_lengths = cp.empty_like(lengths)306 input_batch = batch307 input_lengths = lengths308 309 threads_per_block = 128310 blocks = (n + threads_per_block - 1) // threads_per_block311 312 while next_token < vocab_size and (max_merges is None or merge_count < max_merges):313 # Early stop if all sequences have collapsed (checked on GPU).314 if bool(cp.all(input_lengths == 1)):315 pbar.write("All sequences have collapsed; stopping early.")316 break317 318 factor = next_token # by construction, every token is < next_token319 best = count_best_pair_gpu(input_batch, input_lengths, factor, pad_val)320 if best is None:321 pbar.write("No mergeable pairs found; stopping early.")322 break323 324 best_pair = (best[0], best[1])325 best_freq = best[2]326 if best_freq < 2:327 pbar.write("Best pair frequency is less than 2; stopping early.")328 break329 330 codes.append((best_pair, next_token))331 332 # Launch the merge kernel.333 merge_kernel((blocks,), (threads_per_block,),334 (input_batch,335 work_batch,336 input_lengths,337 work_lengths,338 cp.int64(n),339 cp.int64(L),340 cp.int64(best_pair[0]),341 cp.int64(best_pair[1]),342 cp.int64(next_token),343 cp.int64(pad_val)))344 # Swap buffers for double-buffering.345 input_batch, work_batch = work_batch, input_batch346 input_lengths, work_lengths = work_lengths, input_lengths347 348 next_token += 1349 merge_count += 1350 pbar.update(1)351 pbar.close()352 353 final_batch = cp.asnumpy(input_batch)354 final_lengths = cp.asnumpy(input_lengths)355 final_data = [final_batch[i, :final_lengths[i]].tolist() for i in range(n)]356 return codes, final_data357 358###################################################################################359 360fused_merge_kernel_code = r'''361extern "C" __global__362void fused_merge_kernel(long* data_in, long* data_out, long* lengths, const long pad_val,363 const long num_rows, const long max_len, const long num_merges, const long* merge_rules) {364 int row = blockIdx.x * blockDim.x + threadIdx.x;365 if (row >= num_rows) return;366 long base = row * max_len;367 long cur_len = lengths[row];368 long* cur = data_in + base;369 long* other = data_out + base;370 // Process each merge rule sequentially.371 for (int m = 0; m < num_merges; m++) {372 long a = merge_rules[3 * m];373 long b = merge_rules[3 * m + 1];374 long new_token = merge_rules[3 * m + 2];375 long out_idx = 0;376 for (int i = 0; i < cur_len; i++) {377 if (i < cur_len - 1 && cur[i] == a && cur[i+1] == b) {378 other[out_idx] = new_token;379 out_idx++;380 i++; // Skip the next token.381 } else {382 other[out_idx] = cur[i];383 out_idx++;384 }385 }386 cur_len = out_idx;387 // Swap pointers for the next merge.388 long* temp = cur;389 cur = other;390 other = temp;391 }392 lengths[row] = cur_len;393 // Pad the remaining positions with pad_val.394 for (int i = cur_len; i < max_len; i++) {395 cur[i] = pad_val;396 }397 // If the final result is not in data_in, copy back.398 if (cur != data_in + base) {399 for (int i = 0; i < cur_len; i++) {400 data_in[base + i] = cur[i];401 }402 }403}404'''405 406try:407 fused_kernel = cp.RawKernel(fused_merge_kernel_code, 'fused_merge_kernel')408 409except:410 pass411 412###################################################################################413 414def retokenize_train_data_fused_gpu(train_data, codes, pad_val=-1):415 """416 Retokenize training data using the fully fused GPU kernel.417 418 The entire training dataset is first packed into GPU memory (using pack_sequences).419 All learned merge rules (provided in 'codes') are applied via a single kernel launch.420 Each GPU thread processes one sequence by applying all merge rules sequentially.421 422 Returns:423 tokenized_data: list of retokenized sequences.424 """425 # Pack the data.426 batch, lengths = pack_sequences(train_data, pad_val)427 n, max_len = batch.shape428 # Build a flattened merge_rules array using CuPy.429 if len(codes) > 0:430 merge_rules_list = [[rule[0][0], rule[0][1], rule[1]] for rule in codes]431 merge_rules_gpu = cp.array(merge_rules_list, dtype=cp.int64)432 merge_rules_gpu = merge_rules_gpu.reshape(-1)433 else:434 merge_rules_gpu = cp.empty((0,), dtype=cp.int64)435 num_merges = merge_rules_gpu.shape[0] // 3436 # Preallocate a scratch buffer.437 scratch = cp.empty_like(batch)438 threads_per_block = 128439 blocks = (n + threads_per_block - 1) // threads_per_block440 # Launch the fused kernel.441 fused_kernel((blocks,), (threads_per_block,),442 (batch, scratch, lengths, cp.int64(pad_val),443 cp.int64(n), cp.int64(max_len), cp.int64(num_merges), merge_rules_gpu))444 final_batch = cp.asnumpy(batch)445 final_lengths = cp.asnumpy(lengths)446 tokenized_data = [final_batch[i, :final_lengths[i]].tolist() for i in range(n)]447 return tokenized_data448 449###################################################################################450 451def bpe_encode(seq, codes):452 """453 Iteratively encodes a sequence using BPE merge rules provided in a dictionary.454 455 Args:456 seq (list): A list of tokens (e.g. integers) representing the input sequence.457 codes (dict): A dictionary mapping token pairs (a tuple of two tokens) 458 to a merged token. For example:459 { (1, 2): 100, (100, 3): 101 }460 461 Returns:462 list: The encoded sequence after applying all possible merges.463 464 The function repeatedly scans the entire sequence from left to right;465 whenever it finds a contiguous token pair that exists as a key in the466 codes dict, it replaces that pair with the merged token. This pass is467 repeated until no more merges are possible.468 """469 470 if type(codes) == list:471 codes = dict(codes)472 473 encoded_seq = seq.copy() # work on a copy so as not to modify the original474 done = False475 while not done:476 new_seq = []477 i = 0478 changed = False479 while i < len(encoded_seq):480 # If a merge is possible, merge the two tokens.481 if i < len(encoded_seq) - 1 and (encoded_seq[i], encoded_seq[i + 1]) in codes:482 new_seq.append(codes[(encoded_seq[i], encoded_seq[i + 1])])483 i += 2 # Skip the next token as it was merged.484 changed = True485 else:486 new_seq.append(encoded_seq[i])487 i += 1488 # If no merges occurred in this pass, exit the loop.489 if not changed:490 done = True491 encoded_seq = new_seq492 return encoded_seq493 494###################################################################################495 496def bpe_decode(seq, codes):497 """498 Decodes a sequence encoded with BPE merge rules defined in a codes dictionary.499 500 Args:501 seq (list): The encoded sequence (a list of tokens).502 codes (dict): A dictionary mapping token pairs to the merged token, used during encoding.503 504 Returns:505 list: The fully decoded sequence, with all merged tokens recursively expanded.506 507 The function constructs a reverse mapping that converts a merged token back into 508 its constituent pair. Each token in the sequence is then recursively expanded.509 """510 511 if type(codes) == list:512 codes = dict(codes)513 514 # Build the reverse mapping: key = merged token, value = tuple (original token pair)515 reverse_mapping = {merged: pair for pair, merged in codes.items()}516 517 def recursive_expand(token):518 # If the token is a merged token, expand it recursively.519 if token in reverse_mapping:520 a, b = reverse_mapping[token]521 return recursive_expand(a) + recursive_expand(b)522 else:523 return [token]524 525 decoded_seq = []526 for token in seq:527 decoded_seq.extend(recursive_expand(token))528 return decoded_seq529 530###################################################################################531 532def ensure_triplet(val: Any, name: str = "") -> Tuple[float, float, float]:533 """534 Ensure the given parameter is returned as a triplet.535 If provided as a scalar, promote it to a triplet.536 """537 if np.isscalar(val):538 return (float(val), float(val), float(val))539 elif isinstance(val, (list, tuple)) and len(val) == 3:540 return tuple(float(x) for x in val)541 else:542 raise ValueError(f"{name} must be a scalar or a sequence of 3 numbers.")543 544###################################################################################545 546REP_PENALTY = ensure_triplet(REPETITION_PENALTY, "REPETITION_PENALTY")547SPIKE_STRENGTH = ensure_triplet(SPIKE_PENALTY_STRENGTH, "SPIKE_PENALTY_STRENGTH")548SPIKE_SIG = ensure_triplet(SPIKE_SIGMA, "SPIKE_SIGMA")549 550###################################################################################551 552def sliding_window_view_alternative(a: np.ndarray, window_length: int) -> np.ndarray:553 """554 Create a sliding-window view (without copying) of an array.555 Expected input shape: (n, L, d) and returns: (n, L - window_length + 1, window_length, d)556 """557 n, L, d = a.shape558 new_shape = (n, L - window_length + 1, window_length, d)559 new_strides = (a.strides[0], a.strides[1], a.strides[1], a.strides[2])560 return np.lib.stride_tricks.as_strided(a, shape=new_shape, strides=new_strides)561 562###################################################################################563 564def build_ngram_mapping(data: np.ndarray, memory_len: int) -> Dict[Any, Dict[Any, int]]:565 """566 Build an n-gram mapping from a context (a sequence of triplets) to candidate triplets with frequencies.567 """568 n, L, d = data.shape569 window_length = memory_len + 1 # context (memory) + candidate570 windows = sliding_window_view_alternative(data, window_length)571 # windows shape: (n, L - window_length + 1, window_length, d)572 573 # Split windows into context (first memory_len triplets) and candidates (last triplet)574 contexts = windows[:, :, :memory_len, :] # shape: (n, num_windows, memory_len, d)575 candidates = windows[:, :, memory_len, :] # shape: (n, num_windows, d)576 577 # Flatten the batch and window dimensions.578 contexts_flat = contexts.reshape(-1, memory_len, d)579 candidates_flat = candidates.reshape(-1, d)580 581 mapping = defaultdict(lambda: defaultdict(int))582 total_windows = contexts_flat.shape[0]583 for context_arr, candidate_arr in tqdm.tqdm(584 zip(contexts_flat, candidates_flat),585 total=total_windows,586 desc="Building n-gram mapping"):587 context_key = tuple(map(tuple, context_arr)) # use a tuple of triplets as the key588 candidate_val = tuple(candidate_arr)589 mapping[context_key][candidate_val] += 1590 591 return {context: dict(candidates) for context, candidates in mapping.items()}592 593###################################################################################594 595def precompute_mapping_lookup(mapping: Dict[Any, Dict[Any, int]]) -> Dict[Any, Tuple[Tuple[Any, ...], np.ndarray]]:596 """597 Converts the mapping into a lookup table: context -> (tuple(candidates), frequencies_array).598 """599 mapping_lookup = {}600 for context, candidate_dict in tqdm.tqdm(mapping.items(), desc="Precomputing lookup"):601 candidates = tuple(candidate_dict.keys())602 frequencies = np.array(list(candidate_dict.values()), dtype=np.float64)603 mapping_lookup[context] = (candidates, frequencies)604 return mapping_lookup605 606###################################################################################607 608def build_training_sequences_set(data: np.ndarray) -> set:609 """610 Build a set of training sequences (each as a tuple of triplets) for uniqueness checking.611 """612 return {tuple(map(tuple, seq)) for seq in data}613 614###################################################################################615 616def generate_sequence_optimized(mapping_lookup: Dict[Any, Tuple[Tuple[Any, ...], np.ndarray]],617 training_set: set,618 memory_len: int,619 sequence_length: int = 24,620 max_attempts: int = 1000) -> Optional[Tuple[Tuple[float, float, float], ...]]:621 """622 Autoregressively generate a new, unique sequence using the precomputed mapping lookup.623 The invariant maintained is: the second element of one triplet is never greater than the first element624 of the following triplet.625 626 Two dynamic adjustments are applied for candidate selection:627 628 1. **Dynamic Repetition Penalty:** 629 For each candidate, count the occurrences of each element in the generated sequence.630 Rather than a fixed penalty, this repetition penalty scales with the ratio631 (current_length / sequence_length). In log-space, it subtracts:632 (current_length / sequence_length) * sum_k(count[k] * log(REP_PENALTY[k])633 2. **Dynamic Spike (Variance) Penalty:** 634 For each candidate, compute the squared difference from the running average for each element.635 Use a dynamic sigma that is the maximum between the running standard deviation and the baseline.636 The penalty term for each element is:637 SPIKE_STRENGTH[k] * ((cand[k] - running_avg[k])^2) / (2 * dynamic_sigma[k]^2)638 The overall spike penalty is the sum of the three terms and is subtracted from the candidate’s log frequency.639 640 The resulting candidate log score is computed as:641 log(candidate_frequency) - rep_penalty_component - spike_penalty_component642 A numerical stable softmax is then applied over these scores to determine the probability for drawing a candidate.643 644 If no candidate passing the invariant is found, the attempt is aborted.645 646 Parameters:647 mapping_lookup: Precomputed lookup mapping (context → (candidates, frequencies)).648 training_set: Set of training sequences to ensure uniqueness.649 memory_len: Number of triplets used as context.650 sequence_length: Desired length of the generated sequence.651 max_attempts: Maximum number of generation attempts.652 653 Returns:654 A new unique sequence (tuple of triplets) that respects the invariant, or None if not found.655 """656 mapping_keys = list(mapping_lookup.keys())657 num_keys = len(mapping_keys)658 659 for attempt in range(max_attempts):660 # Select a seed context randomly (from training data so that the invariant holds).661 seed = mapping_keys[np.random.randint(0, num_keys)]662 generated_sequence: List[Tuple[float, float, float]] = list(seed)663 valid_generation = True664 665 while len(generated_sequence) < sequence_length:666 last_triplet = generated_sequence[-1]667 current_context = tuple(generated_sequence[-memory_len:]) # context as tuple of triplets668 candidate_found = False669 670 if current_context in mapping_lookup:671 candidates, frequencies = mapping_lookup[current_context]672 # Filter candidates by invariant:673 # Candidate's first element must be >= last triplet's second element.674 valid_indices = [i for i, cand in enumerate(candidates) if cand[0] >= last_triplet[1]]675 if valid_indices:676 # Filter candidates and their associated frequencies.677 filtered_freqs = frequencies[valid_indices]678 filtered_candidates = [candidates[i] for i in valid_indices]679 680 # Convert candidates into a NumPy array for vectorized operations.681 candidate_array = np.array(filtered_candidates, dtype=np.float64) # shape: (n_candidates, 3)682 683 # Prepare generation history as array.684 generated_array = np.array(generated_sequence, dtype=np.float64) # shape: (T, 3)685 current_length = generated_array.shape[0]686 687 # Running average and standard deviation for dynamic spike adjustment.688 running_avg = np.mean(generated_array, axis=0) # shape: (3,)689 running_std = np.std(generated_array, axis=0) # shape: (3,)690 # Dynamic sigma: ensure a minimum sigma value.691 dynamic_sigma = np.maximum(running_std, np.array(SPIKE_SIG))692 693 # --- Compute Repetition Penalty ---694 # For each candidate, count the number of occurrences for each element along the corresponding column.695 rep_counts = np.array([696 [np.sum(generated_array[:, k] == candidate_array[i, k]) for k in range(3)]697 for i in range(candidate_array.shape[0])698 ]) # shape: (n_candidates, 3)699 # The repetition penalty in log-space.700 rep_penalty_term = np.sum(rep_counts * np.log(np.array(REP_PENALTY)) *701 (current_length / sequence_length), axis=1) # shape: (n_candidates,)702 703 # --- Compute Spike (Variance) Penalty ---704 # Compute the difference per candidate from the running average.705 diff = candidate_array - running_avg # shape: (n_candidates, 3)706 spike_penalty_term = np.sum(np.array(SPIKE_STRENGTH) * (diff**2) / (2 * (dynamic_sigma**2)),707 axis=1) # shape: (n_candidates,)708 709 # --- Compute Candidate Log-Scores ---710 # Use np.log on frequencies (they are positive by construction).711 log_freq = np.log(filtered_freqs)712 log_scores = log_freq - rep_penalty_term - spike_penalty_term713 714 # --- Softmax in Log-space (stable computation) ---715 max_log = np.max(log_scores)716 exp_scores = np.exp(log_scores - max_log)717 probabilities = exp_scores / np.sum(exp_scores)718 719 # Choose the next candidate using advanced probabilities.720 chosen_idx = np.random.choice(len(filtered_candidates), p=probabilities)721 next_triplet = filtered_candidates[chosen_idx]722 candidate_found = True723 724 if not candidate_found:725 # Abort this generation attempt if no valid candidate is available.726 valid_generation = False727 break728 729 generated_sequence.append(next_triplet)730 731 # Ensure the final sequence meets the invariant and is unique.732 if valid_generation and len(generated_sequence) == sequence_length:733 new_sequence = tuple(generated_sequence)734 invariant_ok = all(a[1] <= b[0] for a, b in zip(new_sequence, new_sequence[1:]))735 if invariant_ok and new_sequence not in training_set:736 return new_sequence737 738 return None739 740###################################################################################741 742def analyze_generated_sequence(sequence: tuple, mapping_lookup: dict, memory_len: int) -> tuple:743 """744 Analyze the generated sequence and return several useful statistics745 as both a dictionary and as a nicely formatted string report.746 747 Statistics Computed:748 - unigram_diversity: Ratio of unique triplets to total triplets.749 - repetition_rate: Fraction of repeated triplets.750 - bigram_diversity: Ratio of unique consecutive pairs to total pairs.751 - max_consecutive_repetitions: Maximum number of identical consecutive triplets.752 - avg_candidate_probability (overfit rate): For the transitions (using a sliding window of size753 MEMORY_LEN as context followed by candidate), the average probability of the chosen candidate754 as per the training mapping.755 756 Additional Analytics:757 - element_stats: For each element (index 0, 1, 2) in a triplet, includes:758 * mean, standard deviation, minimum, maximum, and average consecutive absolute difference.759 - avg_transition_entropy: The average entropy of the candidate distributions (from mapping_lookup)760 for each transition context.761 - context_coverage: The fraction of transitions (based on context of length MEMORY_LEN) that are found 762 in the mapping_lookup.763 764 Parameters:765 sequence: Generated sequence (tuple of triplets).766 mapping_lookup: Precomputed mapping lookup.767 memory_len: The context length used.768 769 Returns:770 A tuple containing:771 (stats_dict, stats_report_string)772 """773 stats = {}774 seq_len = len(sequence)775 776 # --- Basic Statistics ---777 778 # Unigram.779 unique_triplets = len(set(sequence))780 stats["unigram_diversity"] = unique_triplets / seq_len781 stats["repetition_rate"] = 1 - (unique_triplets / seq_len)782 783 # Bigram.784 bigrams = [(sequence[i], sequence[i+1]) for i in range(seq_len - 1)]785 unique_bigrams = len(set(bigrams))786 stats["bigram_diversity"] = unique_bigrams / (seq_len - 1)787 788 # Maximum consecutive repetitions.789 max_consecutive = 1790 current_consecutive = 1791 for i in range(1, seq_len):792 if sequence[i] == sequence[i-1]:793 current_consecutive += 1794 if current_consecutive > max_consecutive:795 max_consecutive = current_consecutive796 else:797 current_consecutive = 1798 stats["max_consecutive_repetitions"] = max_consecutive799 800 # Avg Candidate Probability (Overfit Rate)801 overfit_probs = []802 for i in range(memory_len, seq_len):803 context = tuple(sequence[i - memory_len: i])804 candidate = sequence[i]805 if context in mapping_lookup:806 candidates, frequencies = mapping_lookup[context]807 total_freq = np.sum(frequencies)808 try:809 idx = candidates.index(candidate)810 cand_prob = frequencies[idx] / total_freq811 overfit_probs.append(cand_prob)812 except ValueError:813 pass814 stats["avg_candidate_probability"] = np.mean(overfit_probs) if overfit_probs else None815 816 # --- Additional Analytics ---817 818 # 1. Element-Level Statistics.819 seq_arr = np.array(sequence) # shape: (seq_len, 3)820 element_stats = {}821 for dim in range(seq_arr.shape[1]):822 values = seq_arr[:, dim]823 mean_val = np.mean(values)824 std_val = np.std(values)825 min_val = np.min(values)826 max_val = np.max(values)827 # Calculate average absolute difference between consecutive values:828 diffs = np.abs(np.diff(values))829 avg_diff = np.mean(diffs) if diffs.size > 0 else 0830 element_stats[f"element_{dim}"] = {831 "mean": mean_val,832 "std": std_val,833 "min": min_val,834 "max": max_val,835 "avg_consecutive_diff": avg_diff,836 }837 stats["element_stats"] = element_stats838 839 # 2. Transition Entropy:840 entropies = []841 valid_transitions = 0842 for i in range(memory_len, seq_len):843 context = tuple(sequence[i - memory_len: i])844 if context in mapping_lookup:845 candidates, freqs = mapping_lookup[context]846 total_freq = np.sum(freqs)847 if total_freq > 0:848 probs = freqs / total_freq849 # Add a very small constant to avoid log(0)850 epsilon = 1e-10851 entropy = -np.sum(probs * np.log(probs + epsilon))852 entropies.append(entropy)853 valid_transitions += 1854 stats["avg_transition_entropy"] = np.mean(entropies) if entropies else None855 856 # 3. Context Coverage:857 total_transitions = seq_len - memory_len858 stats["context_coverage"] = (valid_transitions / total_transitions) if total_transitions > 0 else None859 860 # --- Build a Pretty Report String ---861 sep_line = "-" * 60862 lines = []863 lines.append(sep_line)864 lines.append("Sequence Analytics Report:")865 lines.append(sep_line)866 lines.append("Overall Statistics:")867 lines.append(f" Unigram Diversity : {stats['unigram_diversity']:.3f}")868 lines.append(f" Repetition Rate : {stats['repetition_rate']:.3f}")869 lines.append(f" Bigram Diversity : {stats['bigram_diversity']:.3f}")870 lines.append(f" Max Consecutive Repetitions: {stats['max_consecutive_repetitions']}")871 cand_prob = stats["avg_candidate_probability"]872 cand_prob_str = f"{cand_prob:.3f}" if cand_prob is not None else "N/A"873 lines.append(f" Avg Candidate Probability : {cand_prob_str}")874 lines.append("")875 876 lines.append("Element-Level Statistics:")877 for dim in sorted(element_stats.keys()):878 ed = element_stats[dim]879 lines.append(f" {dim.capitalize()}:")880 lines.append(f" Mean : {ed['mean']:.3f}")881 lines.append(f" Std Dev : {ed['std']:.3f}")882 lines.append(f" Min : {ed['min']:.3f}")883 lines.append(f" Max : {ed['max']:.3f}")884 lines.append(f" Avg Consecutive Diff : {ed['avg_consecutive_diff']:.3f}")885 lines.append("")886 887 lines.append("Transition Statistics:")888 avg_entropy = stats["avg_transition_entropy"]889 entropy_str = f"{avg_entropy:.3f}" if avg_entropy is not None else "N/A"890 lines.append(f" Average Transition Entropy: {entropy_str}")891 cc = stats["context_coverage"]892 cc_str = f"{cc:.3f}" if cc is not None else "N/A"893 lines.append(f" Context Coverage : {cc_str}")894 lines.append(sep_line)895 896 stats_report = "\n".join(lines)897 898 # Return both the dictionary and the formatted report string.899 return stats, stats_report900 901###################################################################################902 903def autoregressive_generate(start_seq, mel_tones, trg_array, trg_matches_array, num_new_tokens, chunk_len=5):904 905 # Convert sequences to NumPy arrays.906 current_seq = np.array(start_seq, dtype=int) # Shape: (num_tokens, token_dim)907 trg_array = np.array(trg_array, dtype=int) # Shape: (num_candidates, 2, token_dim)908 start_len = len(start_seq)909 910 midx = start_len-1911 912 # Deque for sliding memory of candidate pairs (immutable tuples).913 recent_candidates = deque(maxlen=5)914 915 while (len(current_seq) - start_len) < num_new_tokens:916 917 midx += 1918 919 # Get the last two tokens as context.920 context = current_seq[-(chunk_len-1):] # Shape: (2, token_dim)921 922 sli = 0923 msize = 0924 925 ctx = context[:, :-1].reshape(1, -1)926 trg_mat_arr = trg_matches_array927 928 while msize < 8:929 930 print('=== Slice', sli)931 932 # Compare context with candidates in trg_array.933 match_mask = np.all(ctx == trg_mat_arr, axis=1)934 match_indices = np.where(match_mask)[0]935 936 msize = match_indices.size937 938 if msize < 8:939 sli += 1940 ctx = context[:, :-1].reshape(1, -1)[:, sli:]941 trg_mat_arr = trg_matches_array[:, :-sli]942 943 if match_indices.size == 0:944 if len(current_seq) > start_len:945 946 #tones_chord = sorted([mel_tones[midx], (mel_tones[midx]+7) % 12])947 tones_chord = sorted([mel_tones[midx]])948 new_tuple = [[mel_tones[midx], TMIDIX.ALL_CHORDS_SORTED.index(tones_chord)]] 949 current_seq = np.concatenate((current_seq, new_tuple), axis=0)950 print('Subbed', midx)951 continue952 953 # From the matching candidates, filter out those whose candidate pair is in recent memory.954 available_candidates = []955 cseen = []956 for idx in match_indices:957 958 if idx not in recent_candidates:959 # Convert candidate pair to an immutable tuple960 candidate_pair = tuple(trg_array[idx].tolist())961 if candidate_pair[-1][0] == mel_tones[midx] and candidate_pair[-1][1] not in cseen:962 available_candidates.append((idx, candidate_pair))963 cseen.append(candidate_pair[-1][1])964 965 # If all candidates have recently been used, backtrack.966 if len(available_candidates) < 3:967 if len(current_seq) >= start_len:968 #tones_chord = sorted([mel_tones[midx], (mel_tones[midx]+7) % 12])969 tones_chord = sorted([mel_tones[midx]])970 new_tuple = [[mel_tones[midx], TMIDIX.ALL_CHORDS_SORTED.index(tones_chord)]] 971 current_seq = np.concatenate((current_seq, new_tuple), axis=0)972 #rev_val = random.choice([-1, -2])973 #current_seq = current_seq[:rev_val]974 #print(midx)975 #midx = len(current_seq)976 #print('Reverted', midx, len(current_seq))977 continue978 979 else:980 print(len(available_candidates)) 981 # Choose one available candidate at random.982 chosen_idx, chosen_pair = available_candidates[np.random.choice(len(available_candidates))]983 new_token = trg_array[chosen_idx][-1] # The second token of the candidate pair.984 985 986 # Append the new token to the sequence.987 current_seq = np.concatenate((current_seq, new_token[None, :]), axis=0)988 989 recent_candidates.append(chosen_idx)990 991 print('Gen seq len', len(current_seq))992 993 return current_seq994 995###################################################################################996 997def minkowski_distance_vector_to_matrix(x: cp.ndarray, X: cp.ndarray, p: float = 3) -> cp.ndarray:998 999 """1000 Computes the Minkowski distance between a 1D CuPy array 'x' and each row of a 2D CuPy array 'X'.1001 1002 Parameters:1003 x (cp.ndarray): A 1D array with shape (n_features,) representing a single vector.1004 X (cp.ndarray): A 2D array with shape (n_samples, n_features) where each row is a vector.1005 p (float): The order of the Minkowski distance.1006 For instance:1007 - p=1 yields the Manhattan distance,1008 - p=2 yields the Euclidean distance,1009 - p=3 yields the Minkowski distance and will use the cube-root implementation,1010 - p=∞ (or cp.inf) gives the Chebyshev distance.1011 1012 Returns:1013 cp.ndarray: A 1D array of length n_samples containing the Minkowski distance between 'x' 1014 and the corresponding row in 'X'.1015 """1016 1017 # Compute the element-wise absolute differences between x and every row in X.1018 # Broadcasting x over the rows of X results in an array of shape (n_samples, n_features).1019 diff = cp.abs(X - x)1020 1021 if p == float('inf') or p == cp.inf:1022 # For the Chebyshev distance, use the maximum absolute difference along the feature axis.1023 distances = cp.max(diff, axis=1)1024 elif p == 3:1025 # Instead of using the generic power operation (sum(diff**3) ** (1/3)),1026 # we use cp.cbrt for cube-root calculation when p is exactly 3.1027 distances = cp.cbrt(cp.sum(diff ** 3, axis=1))1028 else:1029 # For general Minkowski distance with finite p,1030 # compute the p-th power of differences, sum them, then take the p-th root.1031 distances = cp.sum(diff ** p, axis=1) ** (1.0 / p)1032 1033 return distances1034 1035###################################################################################1036 1037def pairwise_minkowski_distance(X: cp.ndarray, p: float = 2) -> cp.ndarray:1038 1039 """1040 Computes pairwise Minkowski distances for a 2D CuPy array.1041 1042 Parameters:1043 X (cp.ndarray): A 2D array of shape (n_samples, n_features), where each row represents a vector.1044 p (float): The order of the Minkowski distance.1045 For example:1046 - p=1 is the Manhattan distance,1047 - p=2 is the Euclidean distance,1048 - p=∞ (e.g., float('inf') or cp.inf) is the Chebyshev distance.1049 1050 Returns:1051 cp.ndarray: A 2D array of shape (n_samples, n_samples) containing the pairwise Minkowski distances.1052 """1053 1054 # Use broadcasting to compute the absolute difference between every pair of vectors.1055 # The result of X[:, None, :] - X[None, :, :] will have shape (n_samples, n_samples, n_features).1056 if p == float('inf') or p == cp.inf:1057 # For the Chebyshev distance, take the maximum absolute difference along the feature axis.1058 return cp.max(cp.abs(X[:, None, :] - X[None, :, :]), axis=-1)1059 else:1060 # Raise the absolute differences to the power p.1061 diff_powered = cp.abs(X[:, None, :] - X[None, :, :]) ** p1062 # Sum over the features for each pair (i, j) and then take the p-th root.1063 distances = cp.sum(diff_powered, axis=-1) ** (1.0 / p)1064 1065 return distances1066 1067###################################################################################1068 1069def pairwise_cosine_similarity(X: cp.ndarray, eps: float = 1e-10) -> cp.ndarray:1070 1071 """1072 Computes the pairwise cosine similarity for a 2D CuPy array.1073 1074 Parameters:1075 X (cp.ndarray): A 2D array of shape (n_samples, n_features) where each row represents a vector.1076 eps (float): A small constant added to the denominator to prevent division by zero.1077 1078 Returns:1079 cp.ndarray: A 2D array of shape (n_samples, n_samples) containing the pairwise cosine similarities.1080 """1081 1082 # Compute the dot product between every pair of rows.1083 # This results in a matrix where element (i, j) is the dot product of X[i] and X[j].1084 dot_product = cp.dot(X, X.T)1085 1086 # Compute the L2 norm (Euclidean norm) for each row vector.1087 norms = cp.linalg.norm(X, axis=1)1088 1089 # Compute the outer product of the norms to form the denominator.1090 # The element (i, j) in this matrix is norms[i] * norms[j].1091 norm_matrix = cp.outer(norms, norms)1092 1093 # Compute the cosine similarity matrix.1094 # Adding a small epsilon (eps) to the denominator prevents division by zero.1095 cosine_similarity = dot_product / (norm_matrix + eps)1096 1097 return cosine_similarity1098 1099###################################################################################1100 1101def cosine_similarities(src_array, trg_array):1102 1103 """1104 Computes cosine similarities between 1D src array and 2D trg array1105 """1106 1107 src_norm = cp.linalg.norm(src_array)1108 1109 trg_norms = cp.linalg.norm(trg_array, axis=1)1110 1111 dot_products = cp.dot(trg_array, src_array)1112 1113 cosine_sims = dot_products / (src_norm * trg_norms + 1e-10)1114 1115 return cosine_sims1116 1117###################################################################################1118 1119def embeddings_topk_cosine_neighbors(embeddings,1120 k=1,1121 row_batch=4096,1122 col_batch=40961123 ):1124 1125 """1126 For each embedding, find the indices and similarities of its top-k neighbors,1127 excluding itself, sorted by descending similarity.1128 1129 Args:1130 embeddings (cp.ndarray): shape (N, D), float32 on GPU.1131 k (int): how many neighbors to return (must be < N).1132 row_batch (int): number of rows to process at once.1133 col_batch (int): number of columns to process at once.1134 1135 Returns:1136 top_idx (cp.ndarray): shape (N, k), int32 indices of nearest neighbors.1137 top_sim (cp.ndarray): shape (N, k), float32 cosine similarities.1138 Each row is sorted descending.1139 """1140 1141 # normalize embeddings to unit length1142 norms = cp.linalg.norm(embeddings, axis=1, keepdims=True)1143 embeddings /= norms1144 1145 N, D = embeddings.shape1146 if not (1 <= k < N):1147 raise ValueError(f"k must satisfy 1 ≤ k < N; got N={N}, k={k}")1148 1149 # placeholders for top-k1150 top_sim = cp.full((N, k), -cp.inf, dtype=cp.float32)1151 top_idx = cp.full((N, k), -1, dtype=cp.int32)1152 1153 for i in tqdm.tqdm(range(0, N, row_batch)):1154 i_end = min(i + row_batch, N)1155 rows = embeddings[i:i_end] # (rb, D)1156 rb = i_end - i1157 1158 # per-block buffers1159 best_s = top_sim[i:i_end] # (rb, k)1160 best_i = top_idx[i:i_end] # (rb, k)1161 row_ids = cp.arange(i, i_end) # (rb,)1162 1163 for j in range(0, N, col_batch):1164 j_end = min(j + col_batch, N)1165 cols = embeddings[j:j_end] # (cb, D)1166 sims = rows.dot(cols.T) # (rb, cb)1167 1168 # mask out self-similarities1169 # find rows whose global index ∈ [j, j_end)1170 mask = (row_ids >= j) & (row_ids < j_end)1171 if mask.any():1172 local_rows = cp.where(mask)[0]1173 local_cols = row_ids[mask] - j1174 sims[local_rows, local_cols] = -cp.inf1175 1176 # get top-k within this block1177 # argpartition to grab k largest in each row1178 part = cp.argpartition(sims, -k, axis=1)[:, -k:] # (rb, k)1179 blk_s = sims[cp.arange(rb)[:, None], part] # (rb, k)1180 blk_i = part + j # (rb, k)1181 1182 # merge with running best1183 cat_s = cp.concatenate([best_s, blk_s], axis=1) # (rb, 2k)1184 cat_i = cp.concatenate([best_i, blk_i], axis=1)1185 1186 # select new top-k from the 2k candidates1187 part2 = cp.argpartition(cat_s, -k, axis=1)[:, -k:]1188 best_s = cat_s[cp.arange(rb)[:, None], part2]1189 best_i = cat_i[cp.arange(rb)[:, None], part2]1190 1191 # write back1192 top_sim[i:i_end] = best_s1193 top_idx[i:i_end] = best_i1194 1195 # final sort per row so sims descend1196 if k > 1:1197 order = cp.argsort(-top_sim, axis=1)1198 top_sim = cp.take_along_axis(top_sim, order, axis=1)1199 top_idx = cp.take_along_axis(top_idx, order, axis=1)1200