Constrained tuning: Difference between revisions

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]]''.