###############################################################################
#
#  Temperament finding library -- definitions
#
###############################################################################

import operator, new, string, re
from math import log, exp, sqrt
from UserList import UserList
from UserDict import UserDict

def log2(f):
  return log(f)/log(2)

primeNumbers = (2, 3, 5, 7, 11, 13, 17, 19)

primes = map(log2, primeNumbers[1:])


floatTolerance = 1e-12

###############################################################################
#
#  Stuff for converting between ratios and vectors
#
###############################################################################

def factorizeRatio(n, d=None):
  """return n/d as a ratio space vector"""
  if d==None:
    n, d = n
  assert n>0 and d>0
  return tuple(map(operator.sub, factorize(n), factorize(d)))

def factorize(n):
  assert n>0
  result = []
  residue = n
  for prime in primeNumbers:
    factors = 0
    while residue%prime==0:
      factors = factors + 1
      residue = residue / prime
    result.append(factors)
  assert residue==1
  return tuple(result)

def getRatio(interval):
  """convert the octave-specific ratio space vector
  to a ratio n/d

  return (n,d)

  """
  n=d=1L
  residue = list(interval)
  for i in range(len(interval)):
    while residue[i]>0:
      n = n * primeNumbers[i]
      residue[i] = residue[i]-1L
    while residue[i]<0:
      d = d * primeNumbers[i]
      residue[i] = residue[i]+1L
  try:
    return int(n), int(d)
  except OverflowError:
    return n, d


###############################################################################
#
#  Tonality diamonds
#
###############################################################################


def OddLimit(n):
  harmonics = []
  for n in range(1,n+1,2):
    harmonics.append(factorize(n)[1:])
  return TonalityDiamond(harmonics)

def PrimeDiamond(dimensions, primeHarmonics=None):
  harmonics = [(0,)*dimensions]
  for prime in range(dimensions):
    harmonics.append( (0,)*prime + (1,) + (0,)*(dimensions-prime-1))
  return TonalityDiamond(harmonics, primeHarmonics)

class TonalityDiamond(UserList):
  """all combination of the consonant harmonics"""

  def __init__(self, harmonics=None, primeHarmonics=None):

    if harmonics==None:
      self.data = []
      self.primes = primeHarmonics or []
    elif isinstance(harmonics, TonalityDiamond):
      # this'll be the copy constructor
      self.data = harmonics[:]
      self.harmonics = harmonics.harmonics
      self.primes = primeHarmonics or harmonics.primes
    else:
      if primeHarmonics:
          self.primes = primeHarmonics
          dimensions = len(primeHarmonics)
      else:
        #work out how many primes are being used
        dimensions = 0
        for harmonic in harmonics:
          for i in range(len(harmonic)):
            if i>dimensions and harmonic[i]!=0:
              dimensions = i
        dimensions = dimensions + 1
        self.primes = primes[:dimensions]

      self.harmonics = []
      for harmonic in harmonics:
        self.harmonics.append(harmonic[:dimensions])

      # calculate the diamond
      # or, at least, half of it
      self.data=[]
      for i in range(1,len(harmonics)):
        for j in range(i):
          interval = tuple( map(operator.sub,
                       self.harmonics[i], self.harmonics[j]))
          if not interval in self.data:
            self.data.append(interval)

  def includes(self, interval):
    if tuple(interval) in self.data:
      return 1
    return tuple(map(operator.neg, interval)) in self.data

  def allWithinOctave(self):
    """
    By default, only intervals unique within the
    tritone are enumerated

    This function still leaves out 1:1
    """
    result = []
    for interval in self.data:
      result.append(interval)
      result.append(tuple(map(operator.neg, interval)))
    return result

  def allConsonantChords(self, nNotes):
    root = (0,)*len(self.data[0])
    results = []
    for chord in combinations(
        nNotes-1, # 1/1 is always present
        self.allWithinOctave(),
        self.pairwiseConsonant):
      results.append([root]+chord)
    return results

  def pairwiseConsonant(self, chord):
    for note1, note2 in combinations(2, chord):
      interval = tuple(map(operator.sub, note1, note2))
      if not self.includes(interval):
        return 0
    return 1

  def secondOrder(self):
    """Add and subtract each interval."""
    if not hasattr(self, 'secondOrderCache'):
      result = {}
      for i in range(len(self.data)):
        interval1 = self.data[i]
        for j in range(i+1):
          for fn in (operator.add, operator.sub):
            result[tuple(map(fn,interval1, self.data[j]))]=1
      self.secondOrderCache = result.keys()
    return self.secondOrderCache



limit5 = OddLimit(5)
limit7 = OddLimit(7)
limit9 = OddLimit(9)
limit11 = OddLimit(11)
limit13 = OddLimit(13)
limit15 = OddLimit(15)
limit17 = OddLimit(17)
limit19 = OddLimit(19)
limit21 = OddLimit(21)

simplestOddLimit = map(OddLimit, primeNumbers)


###############################################################################
#
#  Equal Temperaments
#
###############################################################################

def getLimitedETs(consonances, nETs=20, cutoff=0.5, keepDuplicates=0):
  """Get nETs EqualTemperaments consistent in the given consonance limit

  The cutoff is the largest allowable error in terms of step sizes,
  for a consistent mapping.
  0.5 is the value for strict consistency.

  """
  h = consonances.primes
  consist = []
  nNotes=1
  for each in xrange(10000):
    if not len(consist)<nETs:
      return consist
    pet = PrimeET(nNotes, h)
    if cutoff < 0.5:
      ets = [pet]
    else:
      # try all versions where inconsistency is allowed
      ets = alternativeMappings(pet, cutoff)
    for et in ets:
      if et.getWorstError(consonances) < cutoff:
        et.consonances = consonances
        if keepDuplicates or (hcf(et.basis)==1):
          consist.append(et)
    nNotes=nNotes+1
  raise TemperamentException, "not enough ETs found"


class EqualTemperament:
  def __init__(self, basis, targets):
    self.basis = basis
    try:
      self.primes = targets.primes
      self.consonances = targets
    except AttributeError:
      self.primes = targets
      self.consonances = []

  def getWorstError(self, consonances=None):
    """return the largest distance an interval is from JI,
    in terms of scale steps, where each prime interval is
    weighted according to its size.
    """
    consonances = consonances or self.consonances
    errors = self.getErrors()
    worst = 0
    for interval in consonances:
      error = abs(dotprod(interval, errors))
      worst = max(worst, error)
    return worst

  def getErrors(self):
    errors = []
    for index in range(1,len(self.basis)):
      errors.append(self.basis[index] - self.primes[index-1]*self.basis[0])
    return errors

  def getPORMSWE(self):
    """Return the prime, optimum, RMS, weighted error.

    This is the RMS of the prime intervals where octave stretching
    is allowed, with each prime interval weighted according to its size.
    """
    avgStretches, avgSquares = self.getPrimeStretching()
    return sqrt(1.0 - (avgStretches**2 / avgSquares))

  def getPORMSWEStep(self):
    """Return the stretched step size
       for the prime, optimum, RMS, weighted error.
    """
    return self.getPORMSWEStretch() / self.basis[0]

  def getPORMSWEStretch(self):
    """Return the stretch for the prime, optimum, RMS, weighted error.
    """
    avgStretches, avgSquares = self.getPrimeStretching()
    return avgStretches / avgSquares

  def getPrimeStretching(self):
    """Used by getPORMSWE() and getPORMSWEStretch().
    Not likely to be much use on its own.
    """
    sumStretches = sumSquares = 0.0
    for stretch in self.weightedPrimes():
        sumStretches = sumStretches + stretch
        sumSquares = sumSquares + stretch**2
    return sumStretches/len(self.basis), sumSquares/len(self.basis)

  def getTOPError(self):
    """TOP Error (optimal weighted minimax)"""
    wmin, wmax = self.weightedMinMax()
    return (wmax-wmin)/(wmax+wmin)

  def getTOPStretch(self):
    """Return the optimum TOP stretch"""
    wmin, wmax = self.weightedMinMax()
    return 2/(wmin+wmax)

  def weightedMinMax(self):
    """Used for TOP functions.
    Not likely to be much use on its own
    """
    wPrimes = self.weightedPrimes()
    return min(wPrimes), max(wPrimes)

  def weightedPrimes(self):
    """Used for calculating and optimizing weighted prime errors"""
    result = [1.0]
    for i in range(1,len(self.basis)):
      result.append(self.basis[i]/self.primes[i-1]/self.basis[0])
    return result

  def __add__(self, other):
    if self.primes != other.primes:
      raise TypeError("prime intervals must be the same")
    newbasis = []
    for i in range(len(self.basis)):
      newbasis.append(self.basis[i]+other.basis[i])
    return EqualTemperament(newbasis, self.consonances or self.primes)

  def __sub__(self, other):
    if self.primes != other.primes:
      raise TypeError("prime intervals must be the same")
    newbasis = []
    for i in range(len(self.basis)):
      newbasis.append(self.basis[i]-other.basis[i])
    return EqualTemperament(newbasis, self.consonances or self.primes)

  def __and__(self, other):
    if self.primes != other.primes:
      raise TypeError("prime intervals must be the same")
    result = LinearTemperament(
              map(None, self.basis, other.basis),
              self.consonances or other.consonances or self.primes)
    return result;

  def __xor__(self, other):
    return wedgeProduct(self, other)

  def __str__(self):
    return str(self.basis[0])

  def __repr__(self):
    return str(self.basis)

  def __Wedgie__(self):
    return Wedgie(self.basis)


def BestET(nNotes, consonances, cutoff=1.0):
  candidates = []
  for et in alternativeMappings(PrimeET(nNotes, consonances.primes), cutoff):
     candidates.append((et.getWorstError(consonances), et))
  candidates.sort()
  result = candidates[0][1]
  result.consonances = consonances
  return result

def PrimeET(nNotes, primes):
  """return the EqualTemperament with the right
  number of notes to the (formal) octave
  and the nearest approximation to
  each prime interval

  """
  assert nNotes>0
  basis = [nNotes]
  for prime in primes:
    basis.append(nint(prime*nNotes))
  return EqualTemperament(basis, primes)

def WedgieET(wedgie, primesList=None):
    if len(wedgie.keys()[0])>1:
        wedgie = wedgie.complement()
    if not primesList:
        primesList = primes[:wedgie.maxBasis()]
    return(EqualTemperament(wedgie.invariant(), primesList))


def alternativeMappings(et, cutoff, index=1):
  """get all equal temperaments that might be consistent"""

  if index == len(et.basis):
      return [et] ###
  result = alternativeMappings(et, cutoff, index+1)

  target = et.primes[index-1]*et.basis[0]
  approximation = float(et.basis[index])
  error = approximation-target
  if (1.0-abs(error)) < cutoff:
    basis = et.basis[:]
    if error<0:
      basis[index] = basis[index] + 1
    else:
      basis[index] = basis[index] - 1
    result = result + alternativeMappings(
        EqualTemperament(basis, et.primes), cutoff, index+1)

  return result


###############################################################################
#
#  Linear temperaments constructed from equal temperaments or mappings
#
###############################################################################

def defaultFigureOfDemerit(self, consonances=None):
  error = self.getWorstError(consonances)
  try:
    width = self.lastComplexity
  except AttributeError:
    width = self.lastComplexity = self.getComplexity(consonances)
  m, n = self.primesByMelody[0]
  return width*width*error + (
      #fudge factor to give a unique ordering
      1e-10 * (1.618*m+n))


def Temperament(m, n, target):
  """return a Linear Temperament
     consistent with the given ETs
  """
  try:
    et1, et2 = BestET(m, target).basis, BestET(n, target).basis
  except AttributeError:
    et1, et2 = PrimeET(m, target).basis, PrimeET(n, target).basis
  return LinearTemperament(map(None,et1, et2), target)


class LinearTemperament:
  def __init__(self,
      primesByMelody, targets,
      mapping=None, generator=None):

    try:
      primes = targets.primes
      self.consonances = targets
    except AttributeError:
      primes = targets
      self.consonances = simplestOddLimit[len(primes)]

    if isinstance(primesByMelody[0], type(1)):
      # only one equal mapping supplied
      # so we derive the other one from the mapping
      assert mapping

      octave = primesByMelody[0]
      octaveDivision = mapping[0][0]
      period = octave/octaveDivision

      # find another size of octave
      for n in range(period):
        for i in range(len(primesByMelody)):
          if (primesByMelody[i]*n + mapping[i][1])%period:
            break
        else:
          # correct family
          self.primesByMelody = []

          for i in range(len(primesByMelody)):
            steps = (primesByMelody[i]*n + mapping[i][1])/period
            self.primesByMelody.append((primesByMelody[i]-steps, steps))
          break
      else:
        assert generator == 0
        raise TemperamentException, "alternative octave size not found"
    else:
      self.primesByMelody = list(primesByMelody)

    m, n = self.primesByMelody[0]
    octave = m+n
    if m<n:
      m, n = n, m
      for i in range(len(self.primesByMelody)):
        x, y = self.primesByMelody[i]
        self.primesByMelody[i] = (y, x)

    self.primes = primes

    # find the MOS definers
    if mapping:
      octaveDivision = mapping[0][0]
    else:
      octaveDivision = hcf((m,n))
    equivalenceInterval = octave/octaveDivision
    self.generator = generator or getGenerator(n,octave,octaveDivision)

    if mapping:
      self.mapping = mapping
    else:
      self.setMapping()

    self.stdError = 17.0/1200
    # used for measuring complexity
    # the value is for an average listener

    self.FOD = defaultFigureOfDemerit

    # self.basis = 1.0, float(self.generator)/equivalenceInterval
    # nice idea, but it breaks getSmallestContainingMOS()

  def simplify(self):
    # choose a new generator in the light of an optimization
    periodAdjust, newGen = divmod(self.basis[1], self.basis[0])
    periodAdjust = int(periodAdjust)
    # newbasis = oldbasis - periodadjust*period
    if 2*newGen < self.basis[0]:
      self.basis = self.basis[0], newGen
      for i in range(len(self.mapping)):
        equiv, gen = self.mapping[i]
        self.mapping[i] = equiv+gen*int(periodAdjust), gen
    else:
      #newbasis = period-oldbasis+periodadjust*period
      self.basis = self.basis[0], self.basis[0]-newGen
      for i in range(len(self.mapping)):
        equiv, gen = self.mapping[i]
        self.mapping[i] = equiv+gen*(periodAdjust+1), -gen
    period = (
        self.primesByMelody[0][0] +
        self.primesByMelody[0][1]) /self.mapping[0][0]
    generator = self.generator % period
    self.generator = min(generator, period-generator)

  def uncontort(self):
      """Return a temperament like this one but with no contorsion"""
      contorsion = hcf(self.generatorMapping())

      newET = []
      for a, b in self.primesByMelody:
          newET.append(a+b)
      while (hcf(newET+[contorsion]) != 1) :
          for i in range(len(newET)):
              newET[i] = newET[i] / contorsion

      newMapping = []
      torsion = contorsion
      for per, gen in self.mapping:
          torsion = hcf((torsion, per))
          newMapping.append((per, gen/contorsion))
      if torsion != 1:
          # this is the case where the period mapping has a
          # common factor
          # usually (always?) because the two generating ETs
          # had the same common factor
          for i in range(len(newMapping)):
              newMapping[i] = newMapping[i][0] / torsion, newMapping[i][1]

      return LinearTemperament(newET, self.consonances, newMapping)


  def setMapping(self):
      m, n = self.primesByMelody[0]
      octaveDivision = hcf((m,n))

      mGen = getGenerator(n, m, octaveDivision)
      nGen = self.generator - mGen
      determinant = m*nGen-n*mGen
      assert  abs(determinant) == octaveDivision
      determ = determinant/octaveDivision
      equivConv1, equivConv2 = nGen/determ, -mGen/determ
      genConv1, genConv2 = -n/determinant, m/determinant


      # get primes in terms of equivalence interval and generator
      generatorUsed = 0
      self.mapping = [(octaveDivision, 0)]
      for axis1, axis2 in self.primesByMelody[1:]:
        genMapping = genConv1*axis1 + genConv2*axis2
        if genMapping:
          generatorUsed = 1
        self.mapping.append(
            (equivConv1*axis1 + equivConv2*axis2, genMapping))

      if not generatorUsed:
        raise DegenerateException("Intervals are not orthogonal")

  def getFigureOfDemerit(self, consonances=None):
    return self.FOD(self, consonances)

  def setFOD(self, FOD):
    self.FOD = FOD

  def setConsonanceLimit(self, consonances):
    self.consonances = consonances


  def getSmallestContainingMOS(
        self, consonances=None, mustBeProper=0):
    """Dave Keenan's Algorithm"""

    consonances = consonances or self.consonances
    width = self.getWidth(self.getWidestInterval(consonances))

    ratio = abs(self.basis[1]/self.basis[0])
    i = int(ratio)

    mPrev, m = 0, 1
    for each in range(10000):
      if m > width:
        return m*self.mapping[0][0] ###

      ratio = 1/(ratio-i)
      i = int(ratio)

      if mustBeProper:
        mPrev, m = m, m*i + mPrev

      else:
        mPrev, m = m, mPrev

        for n in range(i):
          m = m + mPrev
          if m > width:
            return m*self.mapping[0][0] ###

    raise TemperamentException, "gave up on finding MOS"


  def generatorMapping(self):
    result = []
    for period, generator in self.mapping[1:]:
      result.append(generator)
    return result

  def invariant(self):
    """something that can be used as a unique key"""
    mapping = tuple(self.generatorMapping())
    if filter(None, mapping)[0] < 0:
        new_mapping = []
        for element in mapping:
            new_mapping.append(-element)
        mapping = tuple(new_mapping)
    return (self.mapping[0][0],) + mapping

  def getWidestInterval(self, consonances=None):
    consonances = consonances or self.consonances
    if isinstance(consonances, TonalityDiamond):
      smallest = largest = consonances.harmonics[0]
      smallestWidth = largestWidth = self.getWidth(smallest)
      generatorMapping = self.generatorMapping()
      for interval in consonances.harmonics[1:]:
        intervalWidth = dotprod(interval, generatorMapping)
        if intervalWidth < smallestWidth:
          smallestWidth = intervalWidth
          smallest = interval
        if intervalWidth > largestWidth:
          largestWidth = intervalWidth
          largest = interval
      return map(operator.sub, largest, smallest) ###

    widest = 0
    for interval in consonances:
      width = self.getWidth(interval)
      if width > widest:
        widest = width
        widestInterval = interval
    return widestInterval

  def getWidth(self, interval):
    return abs(dotprod(interval, self.generatorMapping()))

  def getComplexity(self, consonances=None):
    complexity = self.getWidth(self.getWidestInterval(
        consonances))*self.mapping[0][0]
    self.lastComplexity = complexity # used by __cmp__
    return complexity


  def optimizeMinimaxNew(self, consonances):
    """find the minimax using a binary search variant"""

    # initially, choose generators between
    # an unreasonable negative and
    # an unreasonably large interval
    period = 1.0/self.mapping[0][0]
    leftLimit = -10.0
    rightLimit = 10.0 #should catch all but the weirdest
    self.basis = (period, leftLimit)

    rightGrad = 0
    leftValue = 0.0
    lines = []
    for interval in consonances:
      equivs, generators = self.byMapping(interval)
      if generators:
        optimum = dotprod(self.primes, interval)
        leftError = float(equivs)*period - optimum
        lines.append((generators, leftError))
        if abs(generators)>=rightGrad:
          if abs(generators)==rightGrad:
            if abs(leftError) > leftValue:
              # for extreme cases, the worst error will usually be
              # for one of the most complex intervals
              leftValue = abs(leftError)
              rightValue = abs(generators*rightLimit + leftError)
          else:
            rightGrad = abs(generators)
            leftValue = abs(self.getError(interval))
            rightValue = abs(generators*rightLimit + leftError)
    if lines == []:
      raise DegenerateException("no consonances depend on the generator")
    leftGrad = -rightGrad

    for each in xrange(10000):
      # choose a midpoint where the two lines cross
      midPoint = (leftGrad*leftLimit - rightGrad*rightLimit +
          rightValue - leftValue ) / (
          leftGrad - rightGrad)

      """
      find the coeffecients for the lines
      describing the worst error function
      where the generator has its current value
      """
      midValue = 0.0
      midGradL = midGradR = 0
      for errorGradient, axisError in lines:
        error = errorGradient * midPoint + axisError
        if error<0.0:
          errorGradient = -errorGradient
          error = -error

        if abs(error-midValue) < floatTolerance:
          midGradL = min(midGradL, errorGradient)
          midGradR = max(midGradR, errorGradient)
        elif error>midValue:
          midValue = error
          midGradL = midGradR = errorGradient

      if midValue<floatTolerance:
        # not a temperament at all
        self.basis = (period, midPoint)
        return ###

      if midGradL > 0:
        rightLimit, rightGrad, rightValue = midPoint, midGradL, midValue
      else:
        if midGradR > 0:
          # local minimum found
          self.basis = (period, midPoint)
          return ###
        leftLimit, leftGrad, leftValue = midPoint, midGradR, midValue
    raise TemperamentException("minimax overflow for %s" % self.primesByMelody)

  def optimizeMinimaxSimple(self, consonances):
    """simple algorithm that returns a minimax
    but not the preferred one
    """
    generatorMapping = self.generatorMapping()
    optimizingSet = []
    for interval in (consonances):
      if dotprod(interval, generatorMapping):
        optimizingSet.append(interval)
    if len(optimizingSet)==len(consonances):
      optimizingSet = consonances
      # the algorithm is efficient for a TonalityDiamond
    elif optimizingSet==[]:
      raise DegenerateException("degenerate within consonance limit")
    minimax = 1e150
    for interval in optimizingSet:
      self.makeIntervalExact(interval)
      maxError = self.getWorstError(optimizingSet)
      if maxError < minimax:
        minimax = maxError
        bestBasis = self.basis
    try:
      self.basis = bestBasis
    except NameError:
      # this means the minimax is > 1e150 octaves
      print "minimax overflow for %s" % self.primesByMelody


  def optimizeMinimax(self, consonances=None):
    consonances = consonances or self.consonances
    self.optimizeMinimaxNew(consonances)


  def optimizeMinimaxCheck(self, optimizingSet):
    self.optimizeMinimaxSimple(optimizingSet)
    generator = self.basis[1]
    oldOptimum = self.getWorstError(optimizingSet)

    self.optimizeMinimaxNew(optimizingSet)
    if abs(self.basis[1]-generator)>floatTolerance:
      assert oldOptimum > self.getWorstError(optimizingSet)-floatTolerance, "minimaxes %.6f and %.6f differ for %s" % (self.basis[1], generator, self.primesByMelody)

  def makeIntervalExact(self, interval):
    equivsToInterval, gensToInterval = self.byMapping(interval)
    if not gensToInterval:
      raise TemperamentException("doesn't depend on generator")
    genSize = (dotprod(interval, self.primes)
        -float(equivsToInterval)/self.mapping[0][0])/gensToInterval
    self.basis = (1.0/self.mapping[0][0], genSize)

  def byMapping(self, mapped):
    equivsToInterval, gensToInterval = 0, 0
    if len(mapped)==len(self.mapping):
      mapping = self.mapping
    else:
      mapping = self.mapping[1:]
    for prime, primeMap in map(None, mapped, mapping):
      equivsToInterval = equivsToInterval + prime*primeMap[0]
      gensToInterval = gensToInterval + prime*primeMap[1]
    return (equivsToInterval, gensToInterval)

  def getWorstError(self, consonances=None):
    consonances = consonances or self.consonances

    if isinstance(consonances, TonalityDiamond):
      highestError = lowestError = self.getError(consonances.harmonics[0])
      for interval in consonances.harmonics[1:]:
        error = self.getError(interval)
        lowestError = min(error, lowestError)
        highestError = max(error, highestError)
      return highestError-lowestError ###

    worstError = 0
    for interval in consonances:
      error = abs(self.getError(interval))
      worstError = max(worstError, error)
    return worstError

  def getError(self, interval):
    intervalOnBasis = self.byMapping(interval)
    return (
          dotprod(intervalOnBasis, self.basis)
        - dotprod(interval,self.primes))

  def optimizeRMS(self, consonances=None):
    """actually optimize the sum-squared, which comes to the same thing"""

    optimizingSet = consonances or self.consonances

    period = 1/float(self.mapping[0][0])

    numerator = denominator = 0.0
    for interval in optimizingSet:
      equivs, gens = self.byMapping(interval)
      if gens:
        optimum = dotprod(self.primes, interval)
        gradient = float(gens)
        offset = float(equivs)*period - optimum
        #
        # y = mx + c
        #
        # y**2 = m**2*x**2 + 2*cmx + c^^2
        #
        # d(y**2)/dx = 2*m**2*x + 2cm = 0
        #
        # x = -2cm/2m**2 = -cm/m**2
        #
        numerator = numerator - gradient*offset
        denominator = denominator + gradient**2

    self.basis = (period, numerator/denominator)


  def getRMSError(self, consonances=None):
    consonances = consonances or self.consonances

    total = 0.0
    for interval in consonances:
      total = total + self.getError(interval)**2
    return sqrt(total/len(consonances))

  def optimizePORMSWE(self):
    """Set the prime, optimum, RMS, weighted errors optimum"""
    sx0 = sx1 = sx02 = sx12 = sx01 = 0.0

    primes = [1.0]+self.primes
    for i in range(len(primes)):
      m = self.mapping[i]
      x0 = m[0]/primes[i]
      x1 = m[1]/primes[i]
      sx0 = sx0 + x0
      sx1 = sx1 + x1
      sx02 = sx02 + x0**2
      sx12 = sx12 + x1**2
      sx01 = sx01 + x0*x1

    denom = sx02*sx12 - sx01**2
    self.basis = ((sx0*sx12 - sx1*sx01)/denom,
        (sx1*sx02 - sx0*sx01)/denom)

  def getPRMSWError(self):
    """Get the prime, RMS, weighted error"""
    primes = [1.0]+self.primes
    total = 0.0
    for i in range(len(primes)):
      m = self.mapping[i]
      error = (self.basis[0]*m[0] + self.basis[1]*m[1])/primes[i] - 1
      total = total + error*error
    return sqrt(total/len(primes))

  def optimizeTOP(self):
    """Set the TOP generators

    Requires PuLP and GLPK
    """
    import pulp
    prob = pulp.LpProblem("top", pulp.LpMinimize)

    # set three variables: the generators and the thing to minimize
    # all of them have to be positive
    period = pulp.LpVariable("period", 0, None)
    generator = pulp.LpVariable("generator", 0, None)
    error = pulp.LpVariable("error", 0, None)

    # specify that it's the error we want to minimize
    prob.__iadd__((error, "obj"))
    # uses __iadd__() instead of += for syntax compatibility
    # with Python 1.5.2

    # now set the errors of the temperament as constraints
    primes = [1.0]+self.primes
    for i in range(len(primes)):
      weightedPrime = (period*self.mapping[i][0] +
              generator*self.mapping[i][1])/primes[i] 
      # set two error constraints for an overall absolute error
      prob.__iadd__(error >= weightedPrime - 1)
      prob.__iadd__(error >= 1 - weightedPrime)

      prob.solve(pulp.GLPK(msg=0))

      self.basis = period.varValue, generator.varValue

  def getTOPError(self):
    """Get the TOP (not necessarily optimum) error"""
    primes = [1.0]+self.primes
    max = 0.0
    for i in range(len(primes)):
      m = self.mapping[i]
      error = (self.basis[0]*m[0] + self.basis[1]*m[1])/primes[i] - 1
      if abs(error)>max:
          max = abs(error)
    return max

  def getRMSComplexity(self, consonances=None):
    consonances = consonances or self.consonances

    total = 0.0
    ppo = self.mapping[0][0]
    for interval in consonances:
      total = total + (self.getWidth(interval)*ppo)**2
    return sqrt(total/len(consonances))

  def formatEquivalences(self, consonances=None):
    consonances = consonances or self.consonances
    equivalences = self.getEquivalences(consonances)
    if not equivalences:
      equivalences = self.getEquivalences(consonances, order=2, cap=3)
    isharmonic = self.primes == primes[:len(self.primes)]

    results = []
    for set in equivalences:
      line = []
      for interval in set:
        if isharmonic:
          line.append('%i:%i'%getRatio(interval))
        else:
          line.append('(%s)'% string.join(
              map(str, interval), ' '))
      results.append(string.join(line, ' =~ '))

    return string.join(results, '\n')


  def getEquivalences(self, consonances=None, order=1, cap=0):
    """return all sets of intervals that approximate the same"""
    assert order in (1, 2)
    consonances = consonances or self.consonances
    if order==2:
      newConsonances = consonances.secondOrder()
      if cap:
        newConsonances.sort(intervalCompare(self.primes))
        consonances = newConsonances[:len(consonances)*cap]
      else:
        consonances = newConsonances
    intervals = uniqueWithinTritone(consonances)
    octaveDivision = self.mapping[0][0]
    sets = {}
    for interval in intervals:
      mapping = self.byMapping(interval)
      if sets.has_key(mapping):
        sets[mapping].append(interval)
      else:
        sets[mapping] = [interval]
      if octaveDivision%2==0:
        # check for intervals equivalent to a tritone
        if mapping==(octaveDivision/2, 0):
          sets[mapping].append(tuple(map(operator.sub,
              (1,)+(0,)*len(consonances[0]),
              interval)))
    equivalences = []
    for set in sets.values():
      if len(set)>1:
        equivalences.append(set)
    return equivalences

  def getUnisonVectors(self):

    # find some vectors that work
    H = [1]+self.primes
    hcf = self.mapping[0][0]
    fifthIndex = 1
    while self.mapping[fifthIndex][1]==0:
      fifthIndex = fifthIndex + 1
      #okay, so this won't index the fifth any more
    fifth = self.mapping[fifthIndex]
    denom = -hcf*fifth[1]
    vectors = []
    for index in range(1,len(self.mapping)):
      if index == fifthIndex:
        continue
      m, n = self.mapping[index]
      vector = [0]*len(H)
      vector[0] = fifth[1]*m-fifth[0]*n
      vector[fifthIndex] = hcf*n
      vector[index] = denom
      vectors.append(normalizeInterval(vector, H))

    # now simplify them
    nOthers = len(vectors)-1
    cmpfn = intervalCompare(H)
    while nOthers:
      # loop until broken if there are alternatives
      vectors.sort(cmpfn)
      worst = vectors[nOthers]
      for index in range(nOthers):
        alternative = normalizeInterval(
            map(operator.add, vectors[index],worst))
        if cmpfn(alternative, worst)<0:
          vectors[nOthers]=alternative
          break
        alternative = normalizeInterval(
            map(operator.sub, vectors[index],worst), H)
        if cmpfn(alternative, worst)<0:
          vectors[nOthers]=alternative
          break
      else:
        # can't be improved
        break;
    return vectors


  def isUnique(self, consonances=None):
    return not self.getEquivalences(consonances)

  def getScale(self, *rangeArgs, **kwargs):
    """Return a list of cents (by default) values for a generated scale.
    Un-named arguments are as for range()
    and units is the number of steps to the octave

    """
    units = kwargs.get("units", 1200.0)
    scale = []
    period = self.basis[0]
    generator = self.basis[1]
    for nSteps in apply(range, rangeArgs):
      for repeat in range(self.mapping[0][0]):
        root = (nSteps*generator)%period
        scale.append((period*repeat + root)*units)
    scale.sort()
    return scale

  def __repr__(self, consonances=None):
    # send back useful stuff
    consonances = consonances or self.consonances
    if not hasattr(self, 'basis'):
      self.optimizeMinimax()
    params = {
      'generator' : self.generator,
      'basis' : self.basis,
      'mapping' : self.mapping,
      'melody' : self.primesByMelody,
      'width' : self.getWidth(self.getWidestInterval(consonances)),
      'complexity' : self.getComplexity(consonances),
      'mos' : self.getSmallestContainingMOS(consonances),
      'error' : self.getWorstError(consonances),
    }
    if self.primes == primes[:len(self.primes)]:
      units = 1200
      params['unit'] = 'cent'
    else:
      units = 1000
      params['unit'] = 'milli-equivalence'

    params['generatorCents'] = self.basis[1]*units
    params['errorCents'] = self.getWorstError(consonances)*units

    params['period'] = (
        self.primesByMelody[0][0] +
        self.primesByMelody[0][1]) /self.mapping[0][0]

    description = """
%(generator)i/%(period)i, %(generatorCents).1f %(unit)s generator

basis:
%(basis)s

mapping by period and generator:
%(mapping)s

mapping by steps:
%(melody)s

highest interval width: %(width)i
complexity measure: %(complexity)s  (%(mos)s for smallest MOS)
highest error: %(error).6f  (%(errorCents).3f %(unit)ss)""" % params
    if self.isUnique(consonances):
      return description + "\nunique"
    return description


def getLinearTemperaments(ets,
    consonances=None,
    figureOfDemerit=None,
    highestComplexity = 100,
    worstError = None,
    useRMS=0):
  """Returns a list of temperaments as specified,
  ordered by minimax demerit
  """
  consist = ets[:]
  temperaments = TemperamentList()
  while(len(consist)>1):
    et1, consist = consist[0], consist[1:]
    for et2 in consist:
      try:
        tuning = et1 & et2
      except TemperamentException:
        continue
      if consonances:
        tuning.setConsonanceLimit(consonances)
      if figureOfDemerit:
        tuning.setFOD(figureOfDemerit)
      if (tuning.getComplexity()<=highestComplexity
          and hcf(tuning.generatorMapping())==1):
        if useRMS:
          tuning.optimizeRMS(consonances)
          optimum = tuning.getRMSError(consonances)
        else:
          tuning.optimizeMinimax(consonances)
          optimum = tuning.getWorstError(consonances)
        if worstError==None or optimum<worstError:
          key = tuning.invariant()
          if not temperaments.data.has_key(key): # only use the simplest version of each temperament
            tuning.simplify()
            temperaments.add(tuning)

  return temperaments.retrieve()


###############################################################################
#
#  Wedge products
#
###############################################################################

def WedgableRatio(*args):
  return Wedgie(apply(factorizeRatio, args))


class Wedgie(UserDict):

  def __init__(self, seed={}, firstIndex=0):
    if hasattr(seed, "__Wedgie__"):
      other = seed.__Wedgie__()
      self.data = other.data.copy()
    elif isinstance(seed, type({})) or isinstance(seed, UserDict):
      self.data = {}
      for key, value in seed.items():
        self[key] = self[key] + value
    elif isinstance(seed, type(1)) or isinstance(seed, type(1L)) or isinstance(seed, type(1.0)):
        self.data = {():seed}
    elif isinstance(seed, LinearTemperament):
        # wedgies know about linear temperaments
        # so linear temperaments don't have to know about wedgies
        et1 = []
        et2 = []
        for x, y in seed.primesByMelody:
          et1.append(x)
          et2.append(y)
        result = ~(Wedgie(et1)^Wedgie(et2))
        self.data = result.data
    else:
      # assume it's a sequence
      self.data = {}
      dimensions = len(seed)
      # don't return trailing zero elements
      for dimensions in range(len(seed),0,-1):
        if seed[dimensions-1]:
          for i in range(dimensions):
            self[i+firstIndex,]=seed[i]

  def invariant(self, *args):
    wedgie = self.copy()
    wedgie.simplify()
    return apply(wedgie.flatten, args)

  def octaveEquivalent(self):
    result = Wedgie()
    for key, value in self.items():
      if 0 not in key:
        result[key] = value
    return result


  def withinOctave(self, *args, **kwargs):
    assert len(self.keys()[0])==1

    currentPitch = apply(self.pitch, args, kwargs)
    targetPitch = currentPitch%1.0
    result = self.copy()
    result[0,] = result[0,] + nint(targetPitch - currentPitch)

    return result


  def flatten(self, *args):

    basisSize = len(self.keys()[0])
    if basisSize==0:
      return (self.values())
    maxBasis = self.maxBasis()

    # integrity check
    for key in self.keys():
      assert len(key) == basisSize

    if args:
      basisRange = apply(range, args)
    else:
      basisRange = range(maxBasis + 1)

    if basisSize==0:
      return (self.values())

    result = []
    for basis in combinations(basisSize, basisRange):
      result.append(self[tuple(basis)])
    return tuple(result)


  def simplify(self):
    flat = self.flatten()
    gcd = hcf(flat)
    if gcd!=0:
      if filter(None, flat)[0]<0:
        gcd = -gcd
      for key, value in self.items():
        self[key] = value/gcd  #  integer division

  def linearTemperament(self, chroma=None, metric=None, octaveEquivalent=0):

    wedgie = self.copy()
    wedgie.simplify()

    if metric:
      basisDimension = len(metric)+1
      metricList = metric
    else:
      basisDimension = wedgie.maxBasis()+1
      metricList = primes[:basisDimension-1]

    if chroma==None:
      # use the relevant second-order odd limit for candidates
      limit = simplestOddLimit[basisDimension-1].secondOrder()
      allPrimes = [1.0]+metricList
      decorated = []
      for interval in limit:
          vector = reduceWithinTritone(interval, metricList)
          size = dotprod(allPrimes, vector)
          decorated.append((size, vector))
      decorated.sort() # puts small intervals first
      # the first interval in this list is always the unison
      # so we can skip over it
      for size, vector in decorated[1:]:
        chroma = Wedgie(vector)
        if (wedgie^chroma).octaveEquivalent().torsion():
          break
      else:
        chroma = reduceWithinTritone(
            (0,)*(basisDimension-1)+(1,))
        basisDimension = basisDimension + 1
        metricList = primes[:basisDimension-1]

    elif octaveEquivalent:
      chroma = reduceWithinTritone(chroma, metricList)

    generatorMapping = wedgie.octaveEquivalent().complement(1, basisDimension-1)
    if generatorMapping=={}:
            raise TemperamentException, "no generator mapping"
    octaveDivision = generatorMapping.torsion()

    scaleSteps = (wedgie^chroma).complement().invariant(basisDimension)
    conversion = map(None, scaleSteps, generatorMapping.invariant(basisDimension))
    octave = scaleSteps[0]
    period = octave/octaveDivision

    # find the generator
    for generator in range(period):
      mapping = []
      for steps, generators in conversion:
        periods, test = divmod(steps - generator*generators, period)
        if test: break
        mapping.append((periods, generators))
      else:
        return LinearTemperament(
            scaleSteps,
            metricList,
            mapping,
            generator) ###


    raise TemperamentException, "period not found"

  def __invert__(self):
    """ overload ~ operator for the complement"""
    return self.complement()

  def complement(self, firstIndex=0, lastIndex=None):
    """return the equivalent vector with the opposite basis"""

    maxBasis = lastIndex or self.maxBasis()
    result = Wedgie()

    for key, value in self.items():
      newKey = []
      for element in range(firstIndex, maxBasis+1):
        if not element in key:
          newKey.append(element)
      newKey = tuple(newKey)
      result[newKey] = equivalentValue(key+newKey, value)

    return result

  def torsion(self):
    return hcf(self.values())

  def maxBasis(self):
    """return the highest number used in a basis"""
    return max(reduce(operator.add, self.keys(), ()))

  def minBasis(self):
    """return the highest number used in a basis"""
    return min(reduce(operator.add, self.keys(), ()))

  def __getitem__(self, key):
    if not isinstance(key, type(())):
      raise KeyError, "Wedgie keys must be tuples"
    base, value = wedgeEquivalent(key, 1)
    return self.data.get(base, 0)*value

  def __setitem__(self, key, value):
    if not isinstance(key, type(())):
      raise TypeError, "Wedgie keys must be tuples"
    base, value = wedgeEquivalent(key, value)
    if value:
      self.data[base]=value
    else:
      if self.data.has_key(base):
        del self.data[base]

  def __xor__(self, other):
    """overload ^ for the wedge/exterior product"""
    return wedgeProduct(self, other)

  def __rxor__(self, other):
    return wedgeProduct(other, self)

  def __add__(x, b):
    """overload + for elementwise addition"""
    y = makeWedgie(b)
    result = Wedgie()
    for key in x.keys():
      result[key] = result[key] + x[key]
    for key in y.keys():
      result[key] = result[key] + y[key]
    return result

  def __radd__(self, other):
    return makeWedgie(other)+self

  def __sub__(x, b):
    """overload - for elementwise subtraction"""
    y = makeWedgie(b)
    result = Wedgie()
    for element in x.keys()+y.keys():
      result[element]=0
    for key in x.keys():
      result[key] = result[key] + x[key]
    for key in y.keys():
      result[key] = result[key] - y[key]
    return result

  def __rsub__(self, other):
    return makeWedgie(other)-self

  def __neg__(self):
    """overload - for elementwise negation"""
    result = Wedgie()
    for key, value in self.items():
      result[key] = -value
    return result

  def __rmul__(self, other):
    return makeWedgie(other)*self

  def __div__(self, other):
    """overload / for division by a scalar"""
    result = Wedgie()
    for key, value in self.items():
      result[key] = value/other # integer division
    return result

  def __int__(self):
    return int(self[()])

  def __long__(self):
    return long(self[()])

  def __float__(self):
    return float(self[()])

  def __mul__(self, other):
    """overload * for interior product"""
    other = makeWedgie(other)
    try:
      highest = max(self.maxBasis(), other.maxBasis())
      lowest = min(self.minBasis, other.minBasis())
    except ValueError:
      pass
    k = len(self.keys()[0])
    m = len(other.keys()[0])
    n = highest-lowest+1
    product = self^other.complement(lowest, highest)
    result = product.complement(lowest, highest)
    if ((m+k)*(m+k-n))%2:
      return -result
    return result

def wedgeProduct(a, b):
  x = makeWedgie(a)
  y = makeWedgie(b)
  result = Wedgie()
  for base1, value1 in x.items():
    for base2, value2 in makeWedgie(y).items():
      value = value1*value2
      for element in base1:
        if element in base2:
          break
      else:
        base = base1+base2
        result[base] = result[base] + value
  return result


def wedgeEquivalent(base, value):
  """convert the base and value into a
     simpler form consistent with the
     rules of exterior products
  """
  workingBase = list(base)
  for i in range(1, len(base)):
    for j in range(i,0,-1):
      if workingBase[j]<workingBase[j-1]:
        workingBase[j-1], workingBase[j] = (
            workingBase[j],
            workingBase[j-1])
        value = -value
  return tuple(workingBase), value

def equivalentValue(base, value):
    """alternative algorithm that avoids the evil bubblesort"""
    for i in range(len(base)-1):
        for j in range(i+1, len(base)):
            if base[j]<base[i]:
                value = -value
    return value



def makeWedgie(thing):
  if isinstance(thing, Wedgie):
    return thing ###
  return Wedgie(thing)


###############################################################################
#
#  Temperaments generated from unison vectors
#
###############################################################################

def temperOut(unisonVectors, chroma=None, octaveEquivalent = 0):

  if octaveEquivalent:
    unisonVectors = map(reduceWithinTritone, unisonVectors)
  if len(unisonVectors)==1:
    wedgie = Wedgie(unisonVectors[0])
  else:
    wedgie = reduce(wedgeProduct, unisonVectors)

  try:
    missingDimensions = wedgie.maxBasis() - len(unisonVectors) + 1
  except ValueError:
    raise TemperamentException, "linearly dependent"
  if missingDimensions == 1:
    assert chroma==None
    return WedgieET(wedgie) ###
  elif missingDimensions==2:
    return wedgie.linearTemperament(chroma, octaveEquivalent=octaveEquivalent)
  raise TemperamentException, "planar temperaments not supported"


def readVectors(input, octaveEquivalent=0):
  """read some unison vectors from a file-like object

    one set of vectors per line
    octaveEquivalent determines whether the 2 column is returned
    result is a list of lists of tuples
  """
  vectorList = []
  for line in input.readlines():
    if ratioMatch.search(line):
      vectors = []
      for ratio in ratioMatch.findall(line):
        vector = factorizeRatio(map(int, ratio))
        vectors.append(vector[octaveEquivalent:])
      vectorList.append(vectors)

  return vectorList


def getCombinations(vectors):
  dimensions = 0
  ratios = []
  if isinstance(vectors, type("")): # 2.x only:  or isinstance(vectors, type(u"")):
    for ratio in ratioMatch.findall(vectors):
      interval = WedgableRatio(map(long, ratio))
      ratios.append(interval)
      dimensions = max(dimensions, interval.maxBasis()+1)
  else:
    for vector in vectors:
      interval = Wedgie(vector)
      ratios.append(interval)
      dimensions = max(dimensions, interval.maxBasis()+1)
  result = []
  for set in combinations(dimensions-2, ratios):
    result.append(set)
  return result


def unisonVectorsFromIntervals(intervals, cutoff, val=None):
    """
    returns all intervals-between-intervals smaller
    than cutoff
    defaults to normal primes
    input is like the output from TonalityDiamond.secondOrder():
    octave equivalent and only unique within the tritone
    """
    results = []
    val = tuple(val or primes[:len(intervals[0])])
    candidates = {}
    for interval in intervals:
        # catch approx tritones
        candidate = map(operator.add, interval, interval)
        size = dotprod(candidate, val) % 1.0
        if size<cutoff or (1-size)<cutoff:
                candidates[reduceWithinTritone(candidate, val)]=1
        for interval1, interval2 in combinations(2, intervals):
            for func in (operator.add, operator.sub):
                candidate = map(func, interval1, interval2)
                size = dotprod(candidate, val) % 1.0
                if size<cutoff or (1-size)<cutoff:
                          candidates[reduceWithinTritone(candidate, val)]=1
        unison = (0,)*(len(intervals[0])+1)
        if candidates.has_key(unison):
            del candidates[unison]
    return candidates.keys()


def linearTemperamentsFromIntervals(
    intervals,
    consonances=None,
    figureOfDemerit=None,
    highestComplexity = 100,
    worstError = 0.01,
    useRMS=0):
  """
  returns a list of temperaments as specified
  generated by approximating intervals-between-octave
  equivalent intervals to be unisons
  ordered by demerit

  """
  return linearTemperamentsFromUnisonVectors(
      unisonVectorsFromIntervals(intervals, worstError*4))


def linearTemperamentsFromUnisonVectors(
    vectors,
    consonances=None,
    figureOfDemerit=None,
    highestComplexity = 100,
    worstError = 0.01,
    useRMS=0):
  """
  returns a list of temperaments as specified
  generated by approximating octave-specific vectors to be unisons
  ordered by demerit

  """
  temperaments = TemperamentList()
  for set in getCombinations(vectors):
    try:
      tuning = temperOut(set)
      if not isinstance(tuning, EqualTemperament):
        if consonances:
          tuning.setConsonanceLimit(consonances)
        if (tuning.getComplexity()<highestComplexity
              and not temperaments.has_key(tuning.invariant())):
          if figureOfDemerit:
            tuning.setFOD(figureOfDemerit)
          if useRMS:
            tuning.optimizeRMS(consonances)
            optimum = tuning.getRMSError(consonances)
          else:
            tuning.optimizeMinimax(consonances)
            optimum = tuning.getWorstError(consonances)
          if optimum<worstError:
            tuning.simplify()
            temperaments.add(tuning)
    except TemperamentException:
      pass
  return temperaments.retrieve()


ratioMatch = re.compile('(\d+)[/:](\d+)')


def getDimensions(primeLimit):
  for dimension in range(1, len(temper.primeNumbers)):
    if primeLimit == primeNumbers[dimension]:
      return dimension+1
  else:
    raise TemperamentException, "prime limit not supported"


def getDivisors(n):
  """returns the divisors of n in descending order"""
  results = [n]
  for each in range(1,n):
    if n%each==0:
      results.insert(1, each)
  return results


###############################################################################
#
#  Utilities
#
###############################################################################

class TemperamentException(Exception):
  pass

class DegenerateException(TemperamentException):
  pass

def normalizeInterval(vector, metric=None):
  """make the interval ascending"""
  denom = hcf(vector)
  for index in range(len(vector)):
    vector[index] = vector[index]/denom
  if metric and dotprod(vector, metric)<0:
    return map(operator.neg, vector)
  return vector

class intervalCompare:
  def __init__(self, metric):
    self.metric = metric

  def __call__(self, a, b):
    return cmp(self.complexity(a),self.complexity(b))

  def complexity(self, interval):
    """lcm"""
    return dotprod(map(abs, interval), self.metric)


def getGenerator(n, d, hcf=1):
  # Uses Carey&Clampitt's magic fourmula
  gen=1
  while (n*gen+hcf)%d:
    gen = gen+1
    if gen>100000:
      raise TemperamentException, "can't find generator"
  return min(gen, d/hcf - gen)

def dotprod(row, column, add=operator.add, mul=operator.mul):
  return reduce(add, map(mul, row, column))

def hcf(numbers):
  result = 1
  nonZeros = filter(None, map(abs, numbers))
  if len(nonZeros)==0:
    return 0
  if len(nonZeros)==1:
    return nonZeros[0]

  nonZeros.sort()
  result = nonZeros[0]
  residue = nonZeros[1:]
  while residue:
    result = euclidianGCD(result, residue[0])
    residue = residue[1:]
  return result


def euclidianGCD(m, n):
  if m==0: return abs(n)
  if n==0: return abs(m)
  curr, next = m, n
  while next:
    curr, next = next, curr%next
  return curr


def nint(x):
  if x>0:
    return int(x+0.5)
  else:
    return int(x-0.5)


def reduceWithinTritone(interval, basis=None):
  if not basis:
    basis = primes[:len(interval)]
  size = dotprod(interval, basis)
  residue = size%1.0
  if residue>0.5:
    octave = nint(size-residue)
    interval = tuple(map(operator.neg, interval))
  else:
    octave = nint(residue-size)
    interval = tuple(interval)

  while (octave+dotprod(interval, basis)<0.0):
    octave = octave + 1

  return (octave,)+interval


def uniqueWithinTritone(consonances, basis=None):
  results = {}
  if basis==None:
    try:
      basis = consonances.primes
    except AttributeError:
      pass
  for interval in consonances:
    results[reduceWithinTritone(interval, basis)]=1
  return results.keys()


def combinations(number, input, condition=operator.truth):
  """
  Return all combinations of number from the input list.
  The condition is a function returning a boolean to decide
  whether it's worth carrying on.
  This can be useful when a large proportion of the results
  are useless, and this can be predicted early on.
  """
  output = []
  if number==1:
    for next in input:
      if condition([next]):
        output.append([next])
  elif number>0:
    for i in range(len(input)-number+1):
      for element in combinations(number-1, input[i+1:], condition):
        next = [input[i]]+element
        if condition(next):
          output.append(next)
  return output


class TemperamentList:
  def __init__(self):
    self.data = {}
  def add(self, temperament):
    self.data[temperament.invariant()] = (
        temperament.getFigureOfDemerit(),
        temperament)
  def retrieve(self):
    decorated = self.data.values()
    decorated.sort()
    undecorated = []
    for _, temperament in decorated:
      undecorated.append(temperament)
    return undecorated

# we don't need no stinkin' tabs
# vim: tabstop=99