Constrained tuning: Difference between revisions

→Computation: update code
m →Computation: here's a stable version of the code
Line 30: Line 30:


== Computation ==
== Computation ==
As a standard optimization problem, numerous algorithms exist to solve it, such as [[Wikipedia: Sequential quadratic programming|sequential quadratic programming]], to name one. [[Flora Canou]]'s [https://github.com/FloraCanou/temperament_evaluator/blob/85b847a41ab0f6273ea43394147b1b8f8a13b8ea/te_optimizer_legacy.py tuning optimizer] is such an implementation in [https://www.python.org Python]. Note: it uses [https://scipy.org/ Scipy].  
As a standard optimization problem, numerous algorithms exist to solve it, such as [[Wikipedia: Sequential quadratic programming|sequential quadratic programming]], to name one. [[Flora Canou]]'s [https://github.com/FloraCanou/temperament_evaluator/blob/v0.26.4/te_optimizer_legacy.py tuning optimizer] is such an implementation in [https://www.python.org Python]. Note: it uses [https://scipy.org/ Scipy].  


{{Databox| Code |
{{Databox| Code |
<syntaxhighlight lang="python">
<syntaxhighlight lang="python">
# © 2020-2023 Flora Canou | Version 0.26.2
# © 2020-2023 Flora Canou | Version 0.26.4
# 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 43: Line 43:


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]
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 SCALAR:
    CENT = 1200
 
class Norm:  
class Norm:  
    """Norm profile for the tuning space."""
     def __init__ (self, wtype = "tenney", wamount = 1, skew = 0, order = 2):
     def __init__ (self, wtype = "tenney", wamount = 1, skew = 0, order = 2):
         self.wtype = wtype
         self.wtype = wtype
Line 54: Line 57:


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


Line 95: Line 99:
     return linalg.norm (gen @ vals - just_tuning_map, ord = order)
     return linalg.norm (gen @ vals - just_tuning_map, ord = order)


def optimizer_main (vals, subgroup = None, norm = Norm (), #"map" is a reserved word
def optimizer_main (vals, subgroup = None, norm = Norm (),  
         cons_monzo_list = None, des_monzo = None, show = True):
         cons_monzo_list = None, des_monzo = None, show = True):
    # NOTE: "map" is a reserved word
    # optimization is preferably done in the unit of octaves, but for precision reasons
     vals, subgroup = __get_subgroup (vals, subgroup)
     vals, subgroup = __get_subgroup (vals, subgroup)


     just_tuning_map = np.log2 (subgroup)*SCALAR
     just_tuning_map = SCALAR.CENT*np.log2 (subgroup)
     vals_x = norm.weightskewed (vals, subgroup)
     vals_x = norm.weightskewed (vals, subgroup)
     just_tuning_map_x = norm.weightskewed (just_tuning_map, subgroup)
     just_tuning_map_x = norm.weightskewed (just_tuning_map, subgroup)
Line 107: Line 113:
         print ("Euclidean tuning without constraints, solved using lstsq. ")
         print ("Euclidean tuning without constraints, solved using lstsq. ")
     else:
     else:
         gen0 = [SCALAR]*vals.shape[0] #initial guess
         gen0 = [SCALAR.CENT]*vals.shape[0] #initial guess
         cons = () if cons_monzo_list is None else {'type': 'eq', 'fun': lambda gen: (gen @ vals - just_tuning_map) @ cons_monzo_list}
         cons = () if cons_monzo_list is None else {
         res = optimize.minimize (__error, gen0, args = (vals_x, just_tuning_map_x, norm.order), method = "SLSQP",
            'type': 'eq',  
            options = {'ftol': 1e-9}, constraints = cons)
            'fun': lambda gen: (gen @ vals - just_tuning_map) @ cons_monzo_list
        }
         res = optimize.minimize (__error, gen0, args = (vals_x, just_tuning_map_x, norm.order),  
            method = "SLSQP", options = {'ftol': 1e-9}, constraints = cons)
         print (res.message)
         print (res.message)
         if res.success:
         if res.success: