Constrained tuning: Difference between revisions

Get rid of the acronym "ettoc"
Computation: update and + todo
Line 35: Line 35:


== Computation ==
== Computation ==
As a standard optimization problem, numerous algorithms exist to solve it, such as {{w|sequential quadratic programming}}, to name one. [[Flora Canou]]'s [https://github.com/FloraCanou/temperament_evaluator/blob/master/te_optimizer_legacy.py tuning optimizer] solves constrained tuning problems in [https://www.python.org Python], using [https://scipy.org/ Scipy].  
As a standard optimization problem, numerous algorithms exist to solve it, such as {{w|sequential quadratic programming}}, to name one. [[Flora Canou]]'s [https://github.com/FloraCanou/temperament_evaluator Temperament Evaluator] solves constrained tuning problems in [https://www.python.org Python], using [https://scipy.org/ Scipy]'s [https://www.cobyqa.com/stable/ COBYQA] algorithm. Here is an abridged version of it:


{{Todo|rework| make an absolutely minimal version of it. }}
{{Databox| Code |
{{Databox| Code |
<syntaxhighlight lang="python">
<syntaxhighlight lang="python">
# © 2020-2025 Flora Canou
# © 2020-2025 Flora Canou
# This work is licensed under the GNU General Public License version 3.
# This work is licensed under the GNU General Public License version 3.
# Version 0.28.1
# Version 0.30.1


import warnings
import warnings
Line 50: Line 51:
PRIME_LIST = [
PRIME_LIST = [
     2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,  
     2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,  
     41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
     41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89]
]


class SCALAR:
class SCALAR:
Line 59: Line 59:
     """Norm profile for the tuning space."""
     """Norm profile for the tuning space."""


     def __init__ (self, wtype = "tenney", wamount = 1, skew = 0, order = 2):
     def __init__ (self, wtype = None, wmode = 1, wstrength = 1, skew = 0, order = 2):
         self.wtype = wtype
         if wtype:
         self.wamount = wamount
            wmode, wstrength = self.__presets (wtype)
        self.wmode = wmode
         self.wstrength = wstrength
         self.skew = skew
         self.skew = skew
         self.order = order
         self.order = order


     def __get_tuning_weight (self, subgroup):
    @staticmethod
         match self.wtype:
     def __presets (wtype):
         match wtype:  
             case "tenney":
             case "tenney":
                 weight_vec = np.reciprocal (np.log2 (subgroup, dtype = float))
                 wmode, wstrength = 1, 1
             case "wilson" | "benedetti":
             case "wilson" | "benedetti":
                 weight_vec = np.reciprocal (np.array (subgroup, dtype = float))
                 wmode, wstrength = 0, 1
             case "equilateral":
             case "equilateral":
                 weight_vec = np.ones (len (subgroup))
                 wmode, wstrength = 0, 0
             case _:
             case _:
                 warnings.warn ("weighter type not supported, using default (\"tenney\")")
                 warnings.warn ("weighter type not supported, using default (\"tenney\")")
                 self.wtype = "tenney"
                 wmode, wstrength = 1, 1
                return self.__get_weight (subgroup)
         return wmode, wstrength
         return np.diag (weight_vec**self.wamount)


     def __get_tuning_skew (self, subgroup):
     def __weight_vec (self, primes):
        """Returns the interval weight vector for a list of formal primes. """
 
        if not isinstance (self.wmode, (int, np.integer)):
            raise TypeError ("non-integer modes not supported. ")
 
        def modal_weighter (primes, m):
            if m == 0:
                return primes
            elif m > 0:
                return modal_weighter (2*np.log2 (primes), m - 1)
            else:
                return modal_weighter (np.exp2 (primes/2), m + 1)
 
        return (modal_weighter (np.asarray (primes), self.wmode)/2)**self.wstrength
 
    def val_weight (self, primes):
        """Returns the val weight matrix for a list of formal primes. """
        return np.diag (1/self.__weight_vec (primes))
 
    def val_skew (self, subgroup):
        """Returns the val skew matrix for a list of formal primes. """
         if self.skew == 0:
         if self.skew == 0:
             return np.eye (len (subgroup))
             return np.eye (len (subgroup))
Line 91: Line 114:
             r*np.ones ((len (subgroup), 1)), axis = 1)
             r*np.ones ((len (subgroup), 1)), axis = 1)


     def tuning_x (self, main, subgroup):
     def val_transform (self, main, subgroup):
         return main @ self.__get_tuning_weight (subgroup) @ self.__get_tuning_skew (subgroup)
         return main @ self.val_weight (subgroup) @ self.val_skew (subgroup)


def __get_subgroup (main, subgroup):
def __get_subgroup (main, subgroup):
Line 99: Line 122:
         subgroup = PRIME_LIST[:main.shape[1]]
         subgroup = PRIME_LIST[:main.shape[1]]
     elif main.shape[1] != len (subgroup):
     elif main.shape[1] != len (subgroup):
         warnings.warn ("dimension does not match. Casting to the smaller dimension. ")
         warnings.warn ("dimensionalities do not match. Casting to the smaller dimensionality. ")
         dim = min (main.shape[1], len (subgroup))
         dim = min (main.shape[1], len (subgroup))
         main = main[:, :dim]
         main = main[:, :dim]
Line 108: Line 131:
         cons_monzo_list = None, des_monzo = None, show = True):
         cons_monzo_list = None, des_monzo = None, show = True):
     # NOTE: "map" is a reserved word
     # NOTE: "map" is a reserved word
     # optimization is preferably done in the unit of octaves, but for precision reasons
     # optimization would ideally be performed in the unit of octaves
    # unfortunately, that often results in insufficient accuracy
    # the cent is a practical choice of unit, and test shows that further scaling
    # doesn't improve accuracy for most main-sequence temperaments
     breeds, subgroup = __get_subgroup (breeds, subgroup)
     breeds, subgroup = __get_subgroup (breeds, subgroup)


     just_tuning_map = SCALAR.CENT*np.log2 (subgroup)
     just_tuning_map = SCALAR.CENT*np.log2 (subgroup)
     breeds_x = norm.tuning_x (breeds, subgroup)
     breeds_x = norm.val_transform (breeds, subgroup)
     just_tuning_map_x = norm.tuning_x (just_tuning_map, subgroup)
     just_tuning_map_x = norm.val_transform (just_tuning_map, subgroup)
     if norm.order == 2 and cons_monzo_list is None: #simply using lstsq for better performance
     if norm.order == 2 and cons_monzo_list is None: #simply using lstsq for better performance
         res = linalg.lstsq (breeds_x.T, just_tuning_map_x)
         res = linalg.lstsq (breeds_x.T, just_tuning_map_x)
Line 122: Line 148:
         gen0 = just_tuning_map[:breeds.shape[0]] #initial guess
         gen0 = just_tuning_map[:breeds.shape[0]] #initial guess
         if cons_monzo_list is None:
         if cons_monzo_list is None:
             cons = ()
             cons_object = ()
         else:
         else:
             cons = optimize.LinearConstraint ((breeds @ cons_monzo_list).T,  
             cons_object = optimize.LinearConstraint ((breeds @ cons_monzo_list).T,  
                 lb = (just_tuning_map @ cons_monzo_list).T,  
                 lb = (just_tuning_map @ cons_monzo_list).T,  
                 ub = (just_tuning_map @ cons_monzo_list).T)
                 ub = (just_tuning_map @ cons_monzo_list).T)
Line 140: Line 166:
         if np.asarray (des_monzo).ndim > 1 and np.asarray (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 @ breeds @ des_monzo) == 0:
         elif (des_tempered_size := gen @ breeds @ des_monzo) == 0:
             raise ZeroDivisionError ("destretch target is in the nullspace. ")
             raise ZeroDivisionError ("destretch target is in the nullspace. ")
         else:
         else:
             gen *= (just_tuning_map @ des_monzo)/tempered_size
             gen *= (just_tuning_map @ des_monzo)/des_tempered_size


     tempered_tuning_map = gen @ breeds
     tempered_tuning_map = gen @ breeds