User:Hkm/Fokker block code
< User:Hkm
This code can be used to find Fokker blocks.
import numpy as np
from fractions import Fraction
from itertools import combinations, product
from math import gcd
from functools import reduce
from temper.lib_temper.temper import cokernel, defactored_hnf
from temper.lib_temper.optimize import lstsq
from temper.lib_temper.subgroup import p_limit
from temper.lib_temper.interval import factors
def prime_factors(n):
factors = []
for d in range(2, int(n**0.5) + 1):
while n % d == 0: factors.append(d); n //= d
return factors + ([n] if n > 1 else [])
def parse_subgroup(subgroup_str: str) -> list[Fraction]:
if not subgroup_str.strip(): return []
return [Fraction(s) for s in subgroup_str.replace('.', ' ').split()]
def parse_commas(comma_str: str, subgroup: list[Fraction] = [Fraction(p) for p in p_limit(97)]) -> list[Fraction]:
comma_str = comma_str.strip()
if not comma_str:
return []
if comma_str.startswith('['):
try:
cleaned_str = comma_str.replace('[', '').replace(']', '').replace(',', ' ').replace('>', ' >')
monzo_parts = cleaned_str.split('>')
monzo_strs = [part.strip() for part in monzo_parts if part.strip()]
commas = []
for s in monzo_strs:
monzo = np.fromstring(s, dtype=int, sep=' ')
if len(monzo) > len(subgroup):
raise ValueError(f"Monzo length {len(monzo)} does not match subgroup length {len(subgroup)}")
elif len(monzo) < len(subgroup):
monzo = monzo[:len(subgroup)]
comma = reduce(lambda x, y: x * y, [subgroup[i]**monzo[i] for i in range(len(monzo))], 1)
commas.append(comma)
return commas
except ValueError as e:
raise ValueError(f"Failed to parse monzo string: {comma_str}. Error: {e}")
else:
return [Fraction(s) for s in comma_str.split(', ')]
def _setup_temperament_analysis(comma_str: str, subgroup_str: str):
"""Prepares temperament data for analysis."""
if not subgroup_str:
parsed_commas = parse_commas(comma_str)
subgroup_str = ".".join([str(f) for f in sorted(set(
f for c in parsed_commas for f in prime_factors(c.denominator) + prime_factors(c.numerator)))])
subgroup = parse_subgroup(subgroup_str)
else:
subgroup = parse_subgroup(subgroup_str)
parsed_commas = parse_commas(comma_str, subgroup)
prime_subgroup = sorted(list(set(f for s in subgroup for f in prime_factors(s.numerator) + prime_factors(s.denominator))))
if not len(parsed_commas):
mapping = np.eye(len(prime_subgroup))
else:
mapping = defactored_hnf(cokernel(np.hstack([factors(c, prime_subgroup) for c in parsed_commas])))
solution, _ = lstsq((mapping, prime_subgroup), weight="tenney")
return parsed_commas, prime_subgroup, mapping, solution, subgroup_str
def fokker_block(tmonzos: list[tuple], temperament_comma_str: str, subgroup_str: str="",
offset: tuple=(), _prime_subgroup=None, _mapping=None, _solution=None):
""" also allows tmonzos to be a comma string """
if _prime_subgroup is None or _mapping is None or _solution is None:
_, prime_subgroup, mapping, solution, subgroup_str_out = _setup_temperament_analysis(temperament_comma_str, subgroup_str)
else:
prime_subgroup, mapping, solution = _prime_subgroup, _mapping, _solution
subgroup_str_out = subgroup_str
# if tmonzos was passed as a comma string, convert it here. find the commas' monzos and apply the mapping
if isinstance(tmonzos, str):
monzos = [factors(comma, prime_subgroup) for comma in parse_commas(tmonzos, parse_subgroup(subgroup_str_out))]
tmonzos = [np.dot(mapping, monzo) for monzo in monzos]
# forbid enfactored comma bases
reduced_tmonzos = []
for tmonzo in tmonzos:
vector = tmonzo.flatten() # tmonzo is a 2D array, flatten it to get the vector
vector_gcd = reduce(gcd, [abs(x) for x in vector if x != 0])
if vector_gcd > 1: print(f"Warning: enfactored comma bases must be expressed as tmonzos.")
reduced_vector = vector // vector_gcd
reduced_tmonzos.append(tuple(reduced_vector))
tmonzos = reduced_tmonzos
# get list of tmonzos inside parallelotope. offset is left bottom corner. also get tunings
basis_vectors = np.array([t[1:] for t in tmonzos])
dim = basis_vectors.shape[1]
tunings = []
if not offset: offset = np.zeros(dim)
else: offset = np.array(offset)
if basis_vectors.shape[0] != dim: raise ValueError("The number of tmonzo vectors must be"
"equal to the dimension of the space after the first coordinate.")
try: basis_inv = np.linalg.inv(basis_vectors)
except np.linalg.LinAlgError: raise ValueError("tmonzos must be linearly independent.")
iter_min = np.floor(np.sum(np.minimum(0, basis_vectors), axis=0) + offset).astype(int)
iter_max = np.ceil (np.sum(np.maximum(0, basis_vectors), axis=0) + offset).astype(int)
lattice_points = []
ranges = [range(min_c, max_c + 1) for min_c, max_c in zip(iter_min, iter_max)]
for p_tuple in product(*ranges):
p = np.array(p_tuple)
alpha = (p - offset) @ basis_inv
if np.all((alpha >= 0) & (alpha < 1)):
lattice_points.append(p)
points_inside = []
if lattice_points:
s = solution.flatten()
if len(s) != dim + 1:
raise ValueError(f"Number of generators {len(s)} does not match tmonzo dimension {dim+1}.")
for p in lattice_points:
val = np.dot(s[1:], p)
c1 = 1 - np.ceil(val / s[0])
points_inside.append(tuple(np.concatenate(([c1], p)).astype(int)))
tunings.append((val + c1 * s[0]) * 1200) # in cents
# get epimorph val
epimorph_val = cokernel(np.array(tmonzos).T)
# sort points inside by size for checking strength
equal_order = []
zipped = sorted(zip(points_inside, tunings), key=lambda x: x[1])
for point, _ in zipped:
equal_order.append(np.dot(epimorph_val, point))
return zipped, tuple(epimorph_val[0]), all(equal_order[i] == i+1 for i in range(len(equal_order)))
def find_strong_blocks(note_range: tuple[int, int], comma_str: str, subgroup_str: str="", offset: tuple=()):
low_count, high_count = note_range
parsed_commas, prime_subgroup, mapping, solution, subgroup_str_out = _setup_temperament_analysis(comma_str, subgroup_str)
dim = len(prime_subgroup) - len(parsed_commas) - 1
if dim <= 0: return []
s = solution.flatten()
candidate_tmonzos = set()
exponent_ranges = [range(high_count+1)] + [range(-high_count, high_count + 1)] * (dim-1)
for p_tuple in product(*exponent_ranges):
p = np.array(p_tuple)
if np.linalg.norm(p) > high_count:
continue
val = np.dot(s[1:], p)
t_0 = int(np.round(-val / s[0]))
tmonzo = np.concatenate(([t_0], p)).astype(int)
if not np.any(tmonzo):
continue
candidate_tmonzos.add(tuple(tmonzo)) # includes tmonzos with common factors
candidate_tmonzos = sorted(list(candidate_tmonzos), key=lambda x: (sum(abs(c) for c in x), x))
strong_blocks = []
for tmonzo_basis in combinations(candidate_tmonzos, dim):
try:
basis_vectors = np.array([t[1:] for t in tmonzo_basis])
if basis_vectors.shape[0] == basis_vectors.shape[1]:
num_notes = round(abs(np.linalg.det(basis_vectors)))
if not (low_count <= num_notes <= high_count):
continue
result, val, is_strong = fokker_block(
list(tmonzo_basis),
comma_str,
subgroup_str_out,
offset,
_prime_subgroup=prime_subgroup,
_mapping=mapping,
_solution=solution
)
if is_strong:
strong_blocks.append({
'tmonzos': tmonzo_basis,
'notes': result,
'val': val,
'note_count': len(result),
})
except ValueError:
continue
return strong_blocks
if __name__ == "__main__":
# print(*fokker_block([(-2, 20)], "225/224, 1029/1024", offset=(0)), sep="\n")
# print(*fokker_block("25/24", "225/224, 1029/1024", offset=(0)), sep="\n")
# print(*fokker_block("81/80, 49/48", "441/440, 540/539"), sep="\n")
# print(*fokker_block([(-9, 6)], "81/80"), sep="\n")
# print(*fokker_block("25/24", "81/80"), sep="\n")
# print(*fokker_block("1029/1024", "[-4 -1 2>", "2.3.7"), sep="\n")
# print(*fokker_block("1029/1024", "[-4 -1 0 2>"), sep="\n")
# print(*find_strong_blocks((15, 16), "441/440, 540/539"), sep="\n\n")
# print(*find_strong_blocks((9, 9), "", "2.3.5"), sep="\n\n")
print(*find_strong_blocks((4, 26), "49/48", offset=(-2,)), sep="\n\n")