Constrained tuning: Difference between revisions

From Xenharmonic Wiki
Jump to navigation Jump to search
Cmloegcmluin (talk | contribs)
add unchanged-intervals along with eigenmonzo, per prior agreement on Xenwiki Work Group
Update code
Line 34: Line 34:
{{Databox| Code |
{{Databox| Code |
<syntaxhighlight lang="python">
<syntaxhighlight lang="python">
# © 2020-2022 Flora Canou | Version 0.22.1
# © 2020-2023 Flora Canou | Version 0.26.1
# This work is licensed under the GNU General Public License version 3.
# This work is licensed under the GNU General Public License version 3.


Line 42: Line 42:
np.set_printoptions (suppress = True, linewidth = 256, precision = 4)
np.set_printoptions (suppress = True, linewidth = 256, precision = 4)


PRIME_LIST = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61]
PRIME_LIST = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89]
SCALAR = 1200 #could be in octave, but for precision reason
SCALAR = 1200 #could be in octave, but for precision reason


def get_subgroup (main, subgroup):
# norm profile for the tuning space
class Norm:
    def __init__ (self, wtype = "tenney", wamount = 1, skew = 0, order = 2):
        self.wtype = wtype
        self.wamount = wamount
        self.skew = skew
        self.order = order
 
    def __get_weight (self, subgroup):
        if self.wtype == "tenney":
            weight_vec = np.reciprocal (np.log2 (np.array (subgroup, dtype = float)))
        elif self.wtype == "wilson" or self.wtype == "benedetti":
            weight_vec = np.reciprocal (np.array (subgroup, dtype = float))
        elif self.wtype == "equilateral":
            weight_vec = np.ones (len (subgroup))
        else:
            warnings.warn ("weighter type not supported, using default (\"tenney\")")
            self.wtype = "tenney"
            return self.__get_weight (subgroup)
        return np.diag (weight_vec**self.wamount)
 
    def __get_skew (self, subgroup):
        if self.skew == 0:
            return np.eye (len (subgroup))
        elif self.order == 2:
            if self.skew != np.inf:
                r = self.skew/(len (subgroup)*self.skew**2 + 1)
                kr = self.skew*r
            else:
                r = 0
                kr = 1/len (subgroup)
        else:
            raise NotImplementedError ("Weil skew only works with Euclidean norm as of now.")
        return np.append (
            np.eye (len (subgroup)) - kr*np.ones ((len (subgroup), len (subgroup))),
            r*np.ones ((len (subgroup), 1)), axis = 1)
 
    def weightskewed (self, main, subgroup):
        return main @ self.__get_weight (subgroup) @ self.__get_skew (subgroup)
 
def __get_subgroup (main, subgroup):
    main = np.asarray (main)
     if subgroup is None:
     if subgroup is None:
         subgroup = PRIME_LIST[:main.shape[1]]
         subgroup = PRIME_LIST[:main.shape[1]]
Line 55: Line 96:
     return main, subgroup
     return main, subgroup


def get_weight (subgroup, wtype = "tenney", wamount = 1):
def __error (gen, vals, jip, order):
    if wtype == "tenney":
     return linalg.norm (gen @ vals - jip, ord = order)
        weight_vec = np.reciprocal (np.log2 (np.array (subgroup, dtype = float)))
    elif wtype == "wilson" or wtype == "benedetti":
        weight_vec = np.reciprocal (np.array (subgroup, dtype = float))
     elif wtype == "equilateral":
        weight_vec = np.ones (len (subgroup))
    elif wtype == "frobenius":
        warnings.warn ("\"frobenius\" is deprecated. Use \"equilateral\" instead. ")
        weight_vec = np.ones (len (subgroup))
    else:
        warnings.warn ("weighter type not supported, using default (\"tenney\")")
        return get_weight (subgroup, wtype = "tenney", wamount = wamount)
    return np.diag (weight_vec**wamount)


def get_skew (subgroup, skew = 0, order = 2):
def optimizer_main (vals, subgroup = None, norm = Norm (), #"map" is a reserved word
    if skew == 0:
        return np.eye (len (subgroup))
    elif skew == np.inf:
        return np.eye (len (subgroup)) - np.ones ((len (subgroup), len (subgroup)))/len(subgroup)
    elif order == 2:
        r = skew/(len (subgroup)*skew**2 + 1)
        return np.append (
            np.eye (len (subgroup)) - skew*r*np.ones ((len (subgroup), len (subgroup))),
            r*np.ones ((len (subgroup), 1)), axis = 1)
    else:
        raise NotImplementedError ("Weil skew only works with Euclidean norm as of now.")
 
def weightskewed (main, subgroup, wtype = "tenney", wamount = 1, skew = 0, order = 2):
    return main @ get_weight (subgroup, wtype, wamount) @ get_skew (subgroup, skew, order)
 
def error (gen, map, jip, order):
    return linalg.norm (gen @ map - jip, ord = order)
 
def optimizer_main (map, subgroup = None, wtype = "tenney", wamount = 1, skew = 0, order = 2,
         cons_monzo_list = None, des_monzo = None, show = True):
         cons_monzo_list = None, des_monzo = None, show = True):
     map, subgroup = get_subgroup (np.array (map), subgroup)
     vals, subgroup = __get_subgroup (vals, subgroup)


     jip = np.log2 (subgroup)*SCALAR
     jip = np.log2 (subgroup)*SCALAR
     map_wx = weightskewed (map, subgroup, wtype, wamount, skew, order)
     vals_wx = norm.weightskewed (vals, subgroup)
     jip_wx = weightskewed (jip, subgroup, wtype, wamount, skew, order)
     jip_wx = norm.weightskewed (jip, subgroup)
     if order == 2 and cons_monzo_list is None: #te with no constraints, simply use lstsq for better performance
     if norm.order == 2 and cons_monzo_list is None: #simply using lstsq for better performance
         res = linalg.lstsq (map_wx.T, jip_wx)
         res = linalg.lstsq (vals_wx.T, jip_wx)
         gen = res[0]
         gen = res[0]
         print ("L2 tuning without constraints, solved using lstsq. ")
         print ("Euclidean tuning without constraints, solved using lstsq. ")
     else:
     else:
         gen0 = [SCALAR]*map.shape[0] #initial guess
         gen0 = [SCALAR]*vals.shape[0] #initial guess
         cons = () if cons_monzo_list is None else {'type': 'eq', 'fun': lambda gen: (gen @ map - jip) @ cons_monzo_list}
         cons = () if cons_monzo_list is None else {'type': 'eq', 'fun': lambda gen: (gen @ vals - jip) @ cons_monzo_list}
         res = optimize.minimize (error, gen0, args = (map_wx, jip_wx, order), method = "SLSQP",
         res = optimize.minimize (__error, gen0, args = (vals_wx, jip_wx, norm.order), method = "SLSQP",
             options = {'ftol': 1e-9}, constraints = cons)
             options = {'ftol': 1e-9}, constraints = cons)
         print (res.message)
         print (res.message)
Line 112: Line 122:


     if not des_monzo is None:
     if not des_monzo is None:
         if np.array (des_monzo).ndim > 1 and np.array (des_monzo).shape[1] != 1:
         if np.asarray (des_monzo).ndim > 1 and np.asarray (des_monzo).shape[1] != 1:
             raise IndexError ("only one destretch target is allowed. ")
             raise IndexError ("only one destretch target is allowed. ")
         elif (tempered_size := gen @ map @ des_monzo) == 0:
         elif (tempered_size := gen @ vals @ des_monzo) == 0:
             raise ZeroDivisionError ("destretch target is in the nullspace. ")
             raise ZeroDivisionError ("destretch target is in the nullspace. ")
         else:
         else:
             gen *= (jip @ des_monzo)/tempered_size
             gen *= (jip @ des_monzo)/tempered_size


     tuning_map = gen @ map
     tuning_map = gen @ vals
     mistuning_map = tuning_map - jip
     mistuning_map = tuning_map - jip


Line 136: Line 146:
Constraints can be added from the parameter <code>cons_monzo_list</code> of the <code>optimizer_main</code> function. For example, to find the CTE tuning for septimal meantone, you type:  
Constraints can be added from the parameter <code>cons_monzo_list</code> of the <code>optimizer_main</code> function. For example, to find the CTE tuning for septimal meantone, you type:  


<pre>
<syntaxhighlight lang="python">
map = np.array ([[1, 0, -4, -13], [0, 1, 4, 10]])
vals = np.array ([[1, 0, -4, -13], [0, 1, 4, 10]])
optimizer_main (map, cons_monzo_list = np.transpose ([1] + [0]*(map.shape[1] - 1)))
optimizer_main (vals, cons_monzo_list = np.transpose ([1] + [0]*(vals.shape[1] - 1)))
</pre>
</syntaxhighlight>


You should get:  
You should get:  
Line 237: Line 247:


== Systematic name ==
== Systematic name ==
In D&D's guide to RTT, the [[Dave Keenan & Douglas Blumeyer's guide to RTT: alternative complexities #Naming|systematic name]] for the CTE tuning scheme is ''[[Dave Keenan %26 Douglas Blumeyer%27s guide to RTT: all-interval tuning schemes #Held-octave minimax-.28E.29S|held-octave minimax-ES]]'', and the systematic name for the CTWE tuning scheme is ''[[Dave Keenan %26 Douglas Blumeyer%27s guide to RTT: tuning fundamentals #Held-intervals|held-octave]] [[Dave_Keenan %26 Douglas Blumeyer%27s guide to RTT: alternative complexities #Tunings used in 7|minimax-E-lils-S]]''.
In D&D's guide to RTT, the [[Dave Keenan & Douglas Blumeyer's guide to RTT: alternative complexities #Naming|systematic name]] for the CTE tuning scheme is ''[[Dave Keenan %26 Douglas Blumeyer%27s guide to RTT: all-interval tuning schemes #Held-octave minimax-.28E.29S|held-octave minimax-ES]]'', and the systematic name for the CTWE tuning scheme is ''[[Dave Keenan %26 Douglas Blumeyer%27s guide to RTT: tuning fundamentals #Held-intervals|held-octave]] [[Dave_Keenan %26 Douglas Blumeyer%27s guide to RTT: alternative complexities #Tunings used in 7|minimax-E-lils-S]]''.



Revision as of 09:56, 30 May 2023

Constrained tunings are tuning optimization techniques using the constraint of some purely tuned intervals (i.e. eigenmonzos, or unchanged-intervals). CTE tuning (constrained Tenney-Euclidean tuning) is the most typical instance. It has a more sophisticated variant, CTWE tuning (constrained Tenney-Weil-Euclidean tuning). These two tunings will be the focus of this article. Otherwise normed tunings can be defined and computed analogously.

All constrained tunings are standard temperament optimization problems. Specifically, as TE tuning can be viewed as a least squares problem, CTE tuning can be viewed as an equality-constrained least squares problem.

The most common subject of constraint is the octave, which is assumed unless specified otherwise. For higher-rank temperaments, it may make sense to add multiple constraints, such as a pure-2.3 constrained tuning. For a rank-r temperament, specifying m eigenmonzos will yield r - m degrees of freedom to be optimized.

Definition

Given a temperament mapping V and the JIP J, we specify a weight and a skew, represented by transformation matrices W and X, respectively, and a p-norm. Suppose the tuning is constrained by the eigenmonzo list MC. The goal is to find the generator list G by

Minimize

[math]\displaystyle{ \displaystyle \lVert GV_{WX} - J_{WX} \rVert_p }[/math]

subject to

[math]\displaystyle{ \displaystyle (GV - J)M_{\rm C} = O }[/math]

where (·)WX denotes the weight-skew transformation, found by

[math]\displaystyle{ \displaystyle \begin{align} V_{WX} &= VWX \\ J_{WX} &= JWX \end{align} }[/math]

The problem is feasible if

  1. rank (MC) ≤ rank (V), and
  2. The subspaces of MC and N (V) are linearly independent.

Computation

As a standard optimization problem, numerous algorithms exist to solve it, such as sequential quadratic programming, to name one. Flora Canou's tuning optimizer is such an implementation in Python. Note: it depends on Scipy.

Code
# © 2020-2023 Flora Canou | Version 0.26.1
# This work is licensed under the GNU General Public License version 3.

import warnings
import numpy as np
from scipy import optimize, linalg
np.set_printoptions (suppress = True, linewidth = 256, precision = 4)

PRIME_LIST = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89]
SCALAR = 1200 #could be in octave, but for precision reason

# norm profile for the tuning space
class Norm: 
    def __init__ (self, wtype = "tenney", wamount = 1, skew = 0, order = 2):
        self.wtype = wtype
        self.wamount = wamount
        self.skew = skew
        self.order = order

    def __get_weight (self, subgroup):
        if self.wtype == "tenney":
            weight_vec = np.reciprocal (np.log2 (np.array (subgroup, dtype = float)))
        elif self.wtype == "wilson" or self.wtype == "benedetti":
            weight_vec = np.reciprocal (np.array (subgroup, dtype = float))
        elif self.wtype == "equilateral":
            weight_vec = np.ones (len (subgroup))
        else:
            warnings.warn ("weighter type not supported, using default (\"tenney\")")
            self.wtype = "tenney"
            return self.__get_weight (subgroup)
        return np.diag (weight_vec**self.wamount)

    def __get_skew (self, subgroup):
        if self.skew == 0:
            return np.eye (len (subgroup))
        elif self.order == 2:
            if self.skew != np.inf:
                r = self.skew/(len (subgroup)*self.skew**2 + 1)
                kr = self.skew*r
            else:
                r = 0
                kr = 1/len (subgroup)
        else:
            raise NotImplementedError ("Weil skew only works with Euclidean norm as of now.")
        return np.append (
            np.eye (len (subgroup)) - kr*np.ones ((len (subgroup), len (subgroup))),
            r*np.ones ((len (subgroup), 1)), axis = 1)

    def weightskewed (self, main, subgroup):
        return main @ self.__get_weight (subgroup) @ self.__get_skew (subgroup)

def __get_subgroup (main, subgroup):
    main = np.asarray (main)
    if subgroup is None:
        subgroup = PRIME_LIST[:main.shape[1]]
    elif main.shape[1] != len (subgroup):
        warnings.warn ("dimension does not match. Casting to the smaller dimension. ")
        dim = min (main.shape[1], len (subgroup))
        main = main[:, :dim]
        subgroup = subgroup[:dim]
    return main, subgroup

def __error (gen, vals, jip, order):
    return linalg.norm (gen @ vals - jip, ord = order)

def optimizer_main (vals, subgroup = None, norm = Norm (), #"map" is a reserved word
        cons_monzo_list = None, des_monzo = None, show = True):
    vals, subgroup = __get_subgroup (vals, subgroup)

    jip = np.log2 (subgroup)*SCALAR
    vals_wx = norm.weightskewed (vals, subgroup)
    jip_wx = norm.weightskewed (jip, subgroup)
    if norm.order == 2 and cons_monzo_list is None: #simply using lstsq for better performance
        res = linalg.lstsq (vals_wx.T, jip_wx)
        gen = res[0]
        print ("Euclidean tuning without constraints, solved using lstsq. ")
    else:
        gen0 = [SCALAR]*vals.shape[0] #initial guess
        cons = () if cons_monzo_list is None else {'type': 'eq', 'fun': lambda gen: (gen @ vals - jip) @ cons_monzo_list}
        res = optimize.minimize (__error, gen0, args = (vals_wx, jip_wx, norm.order), method = "SLSQP",
            options = {'ftol': 1e-9}, constraints = cons)
        print (res.message)
        if res.success:
            gen = res.x
        else:
            raise ValueError ("infeasible optimization problem. ")

    if not des_monzo is None:
        if np.asarray (des_monzo).ndim > 1 and np.asarray (des_monzo).shape[1] != 1:
            raise IndexError ("only one destretch target is allowed. ")
        elif (tempered_size := gen @ vals @ des_monzo) == 0:
            raise ZeroDivisionError ("destretch target is in the nullspace. ")
        else:
            gen *= (jip @ des_monzo)/tempered_size

    tuning_map = gen @ vals
    mistuning_map = tuning_map - jip

    if show:
        print (f"Generators: {gen} (¢)",
            f"Tuning map: {tuning_map} (¢)",
            f"Mistuning map: {mistuning_map} (¢)", sep = "\n")

    return gen, tuning_map, mistuning_map

optimiser_main = optimizer_main

Constraints can be added from the parameter cons_monzo_list of the optimizer_main function. For example, to find the CTE tuning for septimal meantone, you type:

vals = np.array ([[1, 0, -4, -13], [0, 1, 4, 10]])
optimizer_main (vals, cons_monzo_list = np.transpose ([1] + [0]*(vals.shape[1] - 1)))

You should get:

Optimization terminated successfully.
Generators: [1200.     1896.9521] (¢)
Tuning map: [1200.     1896.9521 2787.8085 3369.5214] (¢)
Mistuning map: [ 0.     -5.0029  1.4948  0.6955] (¢)

Analytical solutions exist for Euclidean (L2) tunings, see Constrained tuning/Analytical solution to constrained Euclidean tunings.

For CTE in particular, it can be solved in the method of Lagrange multiplier. The solution is given by

[math]\displaystyle{ \begin{bmatrix} G^{\mathsf T} \\ \Lambda^{\mathsf T} \end{bmatrix} = \begin{bmatrix} V_{WX} V_{WX}^{\mathsf T} & VM \\ (VM)^{\mathsf T} & O \end{bmatrix}^{-1} \begin{bmatrix} V_{WX} J_{WX}^{\mathsf T}\\ (JM)^{\mathsf T} \end{bmatrix} }[/math]

which is almost an analytical solution. Notice we introduced the vector of lagrange multipliers Λ, with length equal to the number of constraints. The lagrange multipliers have no concrete meaning for the resulting tuning, so they can be discarded.

CTE tuning vs CTWE tuning

Consider the fact that TE tuning does not treat divisive ratios as more important than multiplicative ratios – 5/3 and 15/1 are taken as equally important, for example. To address that, a skew on the space may be introduced, resulting in TWE tuning. Constraining the equave to pure on top of TWE gives CTWE tuning.

POTE tuning works as a quick approximation to CTWE. As POTE destretches the equave, it keeps the angle in the tuning space unchanged, and thus sacrifices multiplicative ratios for divisive ratios. On the contrary, CTE sticks to the original design book of TE as its result remains TE optimal.

These tunings can be very different from each other. Take 7-limit meantone as an example. The POTE tuning map is a little bit flatter than quarter-comma meantone, with all the primes tuned flat:

[math]\displaystyle{ \langle \begin{matrix} 1200.000 & 1896.495 & 2785.980 & 3364.949 \end{matrix} ] }[/math]

The CTWE tuning map is a little bit sharper than quarter-comma meantone, with 5 tuned sharp and 3 and 7 flat:

[math]\displaystyle{ \langle \begin{matrix} 1200.000 & 1896.656 & 2786.625 & 3366.562 \end{matrix} ] }[/math]

The CTE tuning map is even sharper, with 3 tuned flat and 5 and 7 sharp.

[math]\displaystyle{ \langle \begin{matrix} 1200.000 & 1896.952 & 2787.809 & 3369.521 \end{matrix} ] }[/math]

Special constraint

The special eigenmonzo WXj, where j is the all-ones monzo, has the effect of removing the weight-skewed tuning bias. This eigenmonzo is actually proportional to the monzo of the extra dimension introduced by the skew. In other words, it forces the extra dimension to be pure, and therefore, the skew will have no effect with this constrained tuning.

It can be regarded as a distinct optimum. In the case of Tenney weighting, it is the TOCTE tuning (Tenney ones constrained Tenney-Euclidean tuning).

For equal temperaments

Since the one degree of freedom of equal temperaments is determined by the constraint, optimization is not involved. It is thus reduced to TOC tuning (Tenney ones constrained tuning). This constrained tuning demonstrates interesting properties.

The step size g can be found by

[math]\displaystyle{ \displaystyle g = 1/\operatorname {mean} (V_{WX}) }[/math]

The edo number n can be found by

[math]\displaystyle{ \displaystyle n = 1/g = \operatorname {mean} (V_{WX}) }[/math]

Unlike TE or TOP, the optimal edo number space in TOC is linear with respect to V. That is, if V = αV1 + βV2, then

[math]\displaystyle{ \displaystyle \begin{align} n &= \operatorname {mean} (VWX) \\ &= \operatorname {mean} ((\alpha V_1 + \beta V_2)WX) \\ &= \operatorname {mean} (\alpha V_1 WX) + \operatorname {mean} (\beta V_2 WX) \\ &= \alpha n_1 + \beta n_2 \end{align} }[/math]

As a result, the relative error space is also linear with respect to V.

For example, the relative errors of 12ettoc5 (12et in 5-limit TOC) is

[math]\displaystyle{ \displaystyle E_\text {r, 12} = \langle \begin{matrix} -1.55\% & -4.42\% & +10.08\% \end{matrix} ] }[/math]

That of 19ettoc5 is

[math]\displaystyle{ \displaystyle E_\text {r, 19} = \langle \begin{matrix} +4.08\% & -4.97\% & -2.19\% \end{matrix} ] }[/math]

As 31 = 12 + 19, the relative errors of 31ettoc5 is

[math]\displaystyle{ \displaystyle \begin{align} E_\text {r, 31} &= E_\text {r, 12} + E_\text {r, 19} \\ &= \langle \begin{matrix} +2.52\% & -9.38\% & +7.88\% \end{matrix} ] \end{align} }[/math]

Systematic name

In D&D's guide to RTT, the systematic name for the CTE tuning scheme is held-octave minimax-ES, and the systematic name for the CTWE tuning scheme is held-octave minimax-E-lils-S.