#!/usr/bin/env python # -*- coding: utf-8 -*- # *************************************************************************** # * Copyright (C) 2012, Paul Lutus * # * * # * This program is free software; you can redistribute it and/or modify * # * it under the terms of the GNU General Public License as published by * # * the Free Software Foundation; either version 2 of the License, or * # * (at your option) any later version. * # * * # * This program is distributed in the hope that it will be useful, * # * but WITHOUT ANY WARRANTY; without even the implied warranty of * # * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * # * GNU General Public License for more details. * # * * # * You should have received a copy of the GNU General Public License * # * along with this program; if not, write to the * # * Free Software Foundation, Inc., * # * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * # *************************************************************************** import re,sys,os import scipy import matplotlib.pyplot as plt def read_file(path): with open(path) as f: return f.read() def write_file(path,data): with open(path,'w') as f: f.write(data) def ntrp(x,xa,xb,ya,yb): return (x-xa) * (yb-ya) / (xb-xa) + ya def interp(v): global mdata index = 0 result = None top = len(mdata) while(index < top): if(index < top-1 and v >= mdata[index][0] and v <= mdata[index+1][0]): result = index break index += 1 if(result != None): return ntrp(v,mdata[result][0],mdata[result+1][0],mdata[result][1],mdata[result+1][1]) else: return False def format_result(s): v = float(s) return '%24.16e , %24.16e' % (v,interp(v)) # return '%24.16e , %.2f' % (v,interp(v)) # pop program name sys.argv.pop(0) # get number of arguments length = len(sys.argv) if(length < 2): print 'Usage: filename containing measurement data' print ' either: single x value -- to produce single y result,' print ' or: start end step -- to generate result table' quit() path = sys.argv.pop(0) sdata = read_file(path).strip() # get number of arguments again length = len(sys.argv) # convert table data to floats mdata = [] for item in re.findall('(?is)[\d.+\-e]+',sdata): try: mdata.append(float(item)) except: None # create pair array from serial form mdata = [mdata[n:n+2] for n in range(0,len(mdata),2)] # parse x and y data into separate arrays xdata = scipy.array([x[0] for x in mdata]) ydata = scipy.array([x[1] for x in mdata]) # create data pair representation xydata = scipy.array(mdata) try: if(length == 1): print '%s' % compute_show_result(sys.argv[0]) elif(length == 3): a = float(sys.argv.pop(0)) b = float(sys.argv.pop(0)) step = float(sys.argv.pop(0)) oa = 0 s = '' while(a < b): oa = a s = format_result(a) print s a += step if(oa < b): ns = format_result(b) if(ns != s): print ns else: print 'Error: wrong number of arguments.' except Exception as e: print 'Error: %s' % e # create an array of regression results # for use in a graph px = [] py = [] mult = 10 for n in range(int(xdata[-1] * mult)): v = n / mult px.append(v) py.append(interp(v)) pcorr = [] for i,x in enumerate(xdata): pcorr.append(interp(x) - ydata[i]) # saved graphics dots per inch dpiv = 72 a = plt.plot(xdata,ydata,'ro') b = plt.plot(px,py,'b') plt.legend([a,b],['data','prediction'],loc=2) plt.grid(color='g') plt.ylabel('Volume Gallons') plt.xlabel('Sensor Height Inches') plt.title('Linear interpolation regression of %d data points' % len(xdata)) plt.savefig('interp_correlation.png',dpi=dpiv)