TOP vs prime-RMS benchmarks and code

I've written and tested a plain-Python algorithm for TOP now. It's
not the optimal (it looks at n*n points where Gene only suggests
n*(n-1)/2) but isn't much more complicated than the prime-RMS version
(17 lines against 12). The TOP is running more than an order of
magnitude slower and I don't think tweaks to the algorithm will change
that. Of course, if you have better code, I can try it. Here are the
results then, for getting the optimal error for the prime mappings of
all equal temperaments with 10 to 99 steps:

5-limit RMS in 6.8ms
5-limit TOP in 65.3ms
7-limit RMS in 8.3ms
7-limit TOP in 132.0ms
11-limit RMS in 9.5ms
11-limit TOP in 232.7ms
13-limit RMS in 10.4ms
13-limit TOP in 369.7ms
17-limit RMS in 11.5ms
17-limit TOP in 556.2ms
19-limit RMS in 13.2ms
19-limit TOP in 811.0ms

The updated code should be at http://x31eq.com/temper.py
now. I'll copy the relevant methods here -- remember that Yahoo
removes the indentation for web viewing.

There will be faster algorithms, but there's always a speed penalty
for that. The only way I can currently do linear temperaments uses a
library that prints out logging info so there wouldn't be a fair
comparison. Nobody likes that method anyway.

Oh yes, I noticed a bug in the PORMSWE stuff. The step and stretch
methods were mixed up, and one was garbage. So if anybody was using
that, update!

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, stretch):
"""TOP Error for a given octave stretch (non-optimal)"""
worst = 0.0
for w in self.weightedPrimes():
w = w*stretch
if abs(1-w)>worst:
worst = abs(1-w)
return worst

def getTOP(self):
"""Return the TOP error and the optimum stretch"""
best, bestStretch = 1e50, 1.0
for prime1 in self.weightedPrimes():
for prime2 in self.weightedPrimes():
stretch = prime1/prime2
error = self.getTOPError(stretch)
if error < best:
best = error
bestStretch = stretch
return best, bestStretch

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

And here, for reference, is the benchmarking code:

import temper, time
limits = 1, 3, 5, 7, 11, 13, 17, 19
for d in range(2,8):
ets = [temper.PrimeET(n, temper.primes[:d]) for n in range(10,100)]
timestamp = time.time()
for n in range(1000):
for et in ets:
et.getPORMSWE()
print "%2i-limit RMS in %5.1fms" % (limits[d], time.time()-timestamp)
timestamp = time.time()
for n in range(1000):
for et in ets:
et.getTOP()
print "%2i-limit TOP in %5.1fms" % (limits[d], time.time()-timestamp)

Graham