Created
July 17, 2014 01:22
-
-
Save giannitedesco/637a936a91982cfc0c10 to your computer and use it in GitHub Desktop.
Take a list of numbers and find a polynomial function that describes them
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/python | |
| # vim: set fileencoding=utf8 : | |
| import numpy as np | |
| def finite_differences(y): | |
| "Use method of finite differences to determine " | |
| "what order of polynomial to look for" | |
| inp = y[:] | |
| ret = 1 | |
| while len(inp) >= 2: | |
| out = [0 for foo in inp[:-1]] | |
| for j in xrange(1, len(inp)): | |
| out[j - 1] = inp[j] - inp[j - 1] | |
| #print inp | |
| #print out | |
| if len({}.fromkeys(out)) == 1: | |
| return ret | |
| inp = out | |
| ret += 1 | |
| return ret | |
| def sup(s): | |
| "Convert an integer in to a superscript representation in UTF-8" | |
| c = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹'] | |
| out = '' | |
| s = int(s) | |
| while s: | |
| d = s % 10 | |
| s /= 10 | |
| out = c[d] + out | |
| return out | |
| def pretty(p): | |
| "Pretty print a numpy poly1d" | |
| out = '' | |
| for i, c in enumerate(p): | |
| def mul(co, b = True): | |
| if 1.0 == co.round(): | |
| return 'x' | |
| else: | |
| if b: | |
| return '(%g*x)'%co | |
| else: | |
| return '%g*x'%co | |
| power = len(p) - i | |
| if i: | |
| if c < 0: | |
| out = out + ' - ' | |
| c = -c | |
| else: | |
| out = out + ' + ' | |
| if power == 0: | |
| out = out + '%s'%c | |
| elif power == 1: | |
| out = out + mul(c, b = False) | |
| else: | |
| out = out + '%s%s'%(mul(c), sup(power)) | |
| return out | |
| def main(): | |
| "Take a list of numbers and find a polynomial that describes them" | |
| from sys import argv | |
| y = np.array(map(np.float64, argv[1:])) | |
| x = np.array(map(np.float64, range(len(y)))) | |
| order = finite_differences(y) | |
| z = np.polyfit(x, y, order) | |
| #print map(lambda x:'%.3f'%x, z) | |
| p = np.poly1d(z) | |
| #print p | |
| print 'Check results are not too inaccurate:' | |
| print map(lambda x:'%.3f'%p(x), x) | |
| print pretty(p) | |
| if __name__ == '__main__': | |
| try: | |
| main() | |
| raise SystemExit, 0 | |
| except: | |
| raise |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment