#!/usr/local/bin/python

import Numeric, LinearAlgebra, temper, re

class UnisonVectorsException (temper.TemperamentException):
  """For errors specific to unison vector calculations"""

def getLinearTemperament(unisonVectors, includesChroma=1):

  if not includesChroma:
      for i in range(len(unisonVectors)+1):
          chroma = Numeric.zeros(len(unisonVectors)+1)
	  chroma[i]=1
          octaveEquivalent = [tuple(chroma)]+list(unisonVectors)
	  if abs(LinearAlgebra.determinant(
	      Numeric.array(octaveEquivalent)))>0.01:
	      break
	  else:
	      print octaveEquivalent
      else:
          raise UnisonVectorsException, "can't find a chroma"
  else:
    octaveEquivalent = unisonVectors
  octaveSpecific = [(1,)+(0,)*len(octaveEquivalent)]
  octaveSpecific += map(makeOctaveSpecific, octaveEquivalent)

  matrix = Numeric.array(octaveSpecific)
  inverted = adjoint(matrix)
  assert inverted[0][0]>=0
  conversion = Numeric.array([x[:2] for x in inverted])

  # check the conjecture wrt octave equivalent matrices
  oeInverse = adjoint(Numeric.array(octaveEquivalent))
  for index in range(len(octaveEquivalent)):
    if oeInverse[index][0]!=inverted[index+1][1]:
      raise UnisonVectorsException, "two mappings do not agree"

  # divide each column by its gcd
  torsion= temper.hcf([x for x, y in conversion])
  periodsPerOctave = temper.hcf([y for x, y in conversion])
  for row in conversion:
    row /= (torsion, periodsPerOctave)
  octave = conversion[0][0]

  for period in getDivisors(octave):

    # 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:
        if period*periodsPerOctave != octave:
          print "custom period required"

        # correct generator, now make it unique
        if filter(None, mul(mapping, (0, 1)))[0]<1:
          otherMapping = []
          for oldPeriod, oldGenerator in mapping:
            otherMapping.append((oldPeriod + oldGenerator, -oldGenerator))
          mapping = otherMapping

        return temper.LinearTemperament(
            [x for x, y in conversion], temper.primes[:len(octaveEquivalent)],
            mapping, generator)

        raise UnisonVectorsException("alternative octave size not found")
  raise UnisonVectorsException("period not found")

def makeOctaveSpecific(vector):
  primes = temper.primes[:len(vector)]
  size = mul(vector,Numeric.array(primes))
  residue = size%1.0
  if residue>0.5:
    residue = residue-1.0
  octave = residue-size
  return (nint(octave),)+vector


testVectors = [
    [(-1, 2), (4, -1)],
    [(0, -3), (8, 1)],
    [(-1, 2), (8, 1)],
    [(-1, 2), (-4, -2)],
    [(-1, 2), (-2, -1)],
    [(0, -3), (-4, -2)],
    [(3, 1), (7, 6)],
    [(3, 1), (2, -3)],
    [(0, -3, 0), (-4, -2, 0), (-2,0,-1)],
    [(0, -3, 0), (0, 2, -2), (-2, 0, -1)],
    [(-1, 2, 0, 0), (-1, 1, 1, 1), (5, 0, 0, -2), (2, 2, -1, 0)],
    [(-1, 0, 2, 0), (-1, 1, 1, 1), (5, 0, 0, -2), (2, 2, -1, 0)],
    [(2, -1, -1, 0), (-1, 1, 1, 1), (5, 0, 0, -2), (2, 2, -1, 0)],
    [(4, -1, 0, 0), (-1, 1, 1, 1), (5, 0, 0, -2), (2, 2, -1, 0)],
    [(8, 1, 0, 0), (-1, 1, 1, 1), (5, 0, 0, -2), (2, 2, -1, 0)],
    [(-2, 2, 0, -1), (-1, 1, 1, 1), (5, 0, 0, -2), (2, 2, -1, 0)],

    [(2, 2, -1, 0, 0),
     (-29, 0, 0, 0, 0),
     (0, -29, 29, 0, 0),
     (0, 29, 0, -29, 0),
     (0, 0, 0, 29, -29)],


    [(2, 2, -1, 0, 0),
     (2, -1, 2, -1, 0),
     (0, -3, 1, 1, 1),
     (-3, 1, -1, -1, 1),
     (-3, 0, 0, 1, -1)],


    [(0, -3), (4, -1)],

    [(-3, 2, -1,  2,  0),
     ( 5, 0,  0, -2,  0),
     (-3, 0,  0,  1, -1),
     (-3,-2,  0,  0,  2),
     (-1,-2,  4,  0,  0),
     ],

    [(-3,  2, -1,  2,  0),
     ( 5,  0,  0, -2,  0),
     (-1,  0,  0,  0, -1),
     (-4,  2,  0,  0,  0),
     ( 2,  1,  2,  0, -1)],

    [(-3,  2, -1,  2,  0),
     ( 5,  0,  0, -2,  0),
     (-3,  0,  0,  1, -1),
     ( 1, -1,  2, -1,  1),
     ( 1,  0,  2,  0, -1),
    ],

]


mul = Numeric.matrixmultiply

def adjoint(matrix):
  return integerize(
      LinearAlgebra.inverse(matrix)*
      abs(LinearAlgebra.determinant(matrix)))

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

def integerize(matrix):
  return Numeric.array([ map(nint, x)  for x in matrix])

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

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

def getCombinations(text, primeLimit, octaveEquivalent=0):
  """only pairs (7-limit) so far"""
  dimensions = getDimensions(primeLimit)
  ratios = [
      tuple(temper.factorizeRatio(*map(int,x))[octaveEquivalent:dimensions])
      for x in ratioMatch.findall(text)]
  result = []
  for i in range(len(ratios)-1):
    for ratio in ratios[i+1:]:
      result.append((ratios[i],ratio))
  return result

def readVectors(input, primeLimit, chromaticLast=0, octaveEquivalent=0):
  """read some unison vectors from a file-like object
  
    the prime limit determines the size of the result

    chromaticLast=1 reverses the order of the vectors.
    Set it if the chromatic UV comes at the end of the list
    (default is at the start
  
    octaveEquivalent determines whether the 2 column is returned
  
    result is a list of lists of tuples
  """
  vectorList = []
  dimensions = getDimensions(primeLimit)
  line = input.readline()
  while 1:
    vectors = []
    while ratioMatch.search(line)==None:
      if line=='': return vectorList ###
      line = input.readline()
    while ratioMatch.search(line):
      for ratio in ratioMatch.findall(line):
        vector = temper.factorizeRatio(*map(int, ratio))
        vectors.append(tuple(vector[octaveEquivalent:dimensions]))
        assert filter(None, vector[dimensions:])==[]
      line = input.readline()
    if chromaticLast:
      vectors.reverse()
    vectorList.append(vectors)

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


if __name__=='__main__':
  for unisonVectors in testVectors:
    vectors = Numeric.array(unisonVectors)
    print "\n\nUnison vectors"
    print vectors
    print "\n\nadjoint of that matrix"
    inv = adjoint(vectors)
    print inv
    try:
      scale = getLinearTemperament(unisonVectors)
      print scale
      firstMapping = scale.mapping
      # see if you can make do without the chromatic unison vector
      nPrimes = len(unisonVectors[0])
      for chromaticPrime in range(nPrimes):
        chroma = [0]*nPrimes
        chroma[chromaticPrime]=1
        unisonVectors[0]=tuple(chroma)
        vectors = Numeric.array(unisonVectors)
        try:
          scale = getLinearTemperament(unisonVectors)
          if scale.mapping != firstMapping:
            print "discrepancy"
            print firstMapping
            print scale.mapping
          break
        except LinearAlgebra.LinAlgError:
          pass
        except UnisonVectorsException, uve:
          print uve
      else:
        print "chroma not found"
    except temper.TemperamentException, te:
      print te