Constrained tuning: Difference between revisions

Versus POTE tuning: +more discussion on the difference
Computation: code update
Line 44: Line 44:
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.  
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.  


Otherwise, 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/b47462bff4a6146e46ee82b4425b1d2410c00e46/te_optimizer_legacy.py tuning optimizer] is such an implementation in [https://www.python.org Python]. Note: it depends on [https://scipy.org/ Scipy].  
Otherwise, 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/55859971e5db037f5353267c328ca4a20f515ecd/te_optimizer_legacy.py tuning optimizer] is such an implementation in [https://www.python.org Python]. Note: it depends on [https://scipy.org/ Scipy].  


{{Databox| Code |
{{Databox| Code |
<syntaxhighlight lang="python">
<syntaxhighlight lang="python">
# © 2020-2022 Flora Canou | Version 0.15.1
# © 2020-2022 Flora Canou | Version 0.20.0
# This work is licensed under the GNU General Public License version 3.
# This work is licensed under the GNU General Public License version 3.


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


Line 59: Line 59:
SCALAR = 1200 #could be in octave, but for precision reason
SCALAR = 1200 #could be in octave, but for precision reason


def subgroup_normalize (main, subgroup):
def get_subgroup (main, subgroup):
     if subgroup is None:
     if subgroup is None:
         subgroup = PRIME_LIST[:main.shape[1]]
         subgroup = PRIME_LIST[:main.shape[1]]
Line 69: Line 69:
     return main, subgroup
     return main, subgroup


def weighted (matrix, subgroup, wtype = "tenney"):
def get_weight (subgroup, wtype = "tenney", *, k = 0.5):
    if not wtype in {"tenney", "frobenius", "partch"}:
        wtype = "tenney"
        warnings.warn ("unknown weighter type, using default (\"tenney\")")
 
     if wtype == "tenney":
     if wtype == "tenney":
         weighter = np.diag (1/np.log2 (subgroup))
         return np.diag (1/np.log2 (subgroup))
     elif wtype == "frobenius":
     elif wtype == "frobenius":
         weighter = np.eye (len (subgroup))
         return np.eye (len (subgroup))
     elif wtype == "partch":
     elif wtype == "inverse tenney":
         weighter = np.diag (np.log2 (subgroup))
        return np.diag (np.log2 (subgroup))
     return matrix @ weighter
    elif wtype == "benedetti":
         return np.diag (1/np.array (subgroup))
    elif wtype == "weil":
        return linalg.pinv (np.append ((1 - k)*np.diag (np.log2 (subgroup)), [k*np.log2 (subgroup)], axis = 0))
     else:
        warnings.warn ("weighter type not supported, using default (\"tenney\")")
        return get_weight (subgroup, wtype = "tenney")
 
def weighted (main, subgroup, wtype = "tenney", *, k = 0.5):
    return main @ get_weight (subgroup, wtype = wtype, k = k)


def error (gen, map, jip, order = 2):
def error (gen, map, jip, order = 2):
     return linalg.norm (gen @ map - jip, ord = order)
     return linalg.norm (gen @ map - jip, ord = order)


def optimizer_main (map, subgroup = None, wtype = "tenney", order = 2, cons_monzo_list = None, stretch_monzo = None, show = True):
def optimizer_main (map, subgroup = None, wtype = "tenney", order = 2,
     map, subgroup = subgroup_normalize (np.array (map), subgroup)
        cons_monzo_list = None, des_monzo = None, show = True, *, k = 0.5):
     map, subgroup = get_subgroup (np.array (map), subgroup)
 
    if wtype == "weil" and order != 2:
        warnings.warn ("The weil weighter as of now is only meant to be used for L2.")


     jip = np.log2 (subgroup)*SCALAR
     jip = np.log2 (subgroup)*SCALAR
     map_w = weighted (map, subgroup, wtype = wtype)
     map_w = weighted (map, subgroup, wtype = wtype, k = k)
     jip_w = weighted (jip, subgroup, wtype = wtype)
     jip_w = weighted (jip, subgroup, wtype = wtype, k = k)
     if order == 2 and cons_monzo_list is None: #te with no constraints, simply use lstsq for better performance
     if order == 2 and cons_monzo_list is None: #te with no constraints, simply use lstsq for better performance
         res = linalg.lstsq (map_w.T, jip_w)
         res = linalg.lstsq (map_w.T, jip_w)
Line 105: Line 114:
             raise ValueError ("infeasible optimization problem. ")
             raise ValueError ("infeasible optimization problem. ")


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


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


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


     return gen, tuning_map
     return gen, tuning_map, mistuning_map


optimiser_main = optimizer_main
optimiser_main = optimizer_main