""" ToyA.py Makes plots from data created by the GUI program Hview.py EG: python ToyA.py Save_Hview_T_100 Save_Hview_T_118 Save_Hview_T_119 \ figs/ToyTS.pdf figs/ToyStretch.pdf Copyright (c) 2005, 2007 Andrew Fraser This file is part of HMM_DS_Code. HMM_DS_Code 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. HMM_DS_Code 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 sys, math, numpy, matplotlib, numpy.linalg EIG = numpy.linalg.eigh INV = numpy.linalg.inv matplotlib.use('PDF') import pylab assert len(sys.argv) == 6,"""%s: Wrong number of arguments. Syntax: python %s dataA dataB dataC fig0 fig1"""%(sys.argv[0],sys.argv[0]) dataA, dataB, dataC, fig0, fig1 = sys.argv[1:6] params = {'axes.labelsize': 12, 'text.fontsize': 10, 'legend.fontsize': 10, 'text.usetex': True, 'xtick.labelsize': 11, 'ytick.labelsize': 11} pylab.rcParams.update(params) def read_TS(name): """ Read a four component time series from the file with name "name". Start with the first sample and continue for 101 samples. The file should have been written by Hview.py """ lines = open(name,'r').readlines() assert lines[1].split()[0] == 'Tmin=' and lines[1].split()[2] == 'Tmax=',""" unexpected format in file %s missed "Tmin=" or "Tmax=" in line 1: %s"""%(name,lines[1]) Tmin = int(lines[1].split()[1]) Tmax = int(lines[1].split()[3]) Tmax = min(Tmax,Tmin+101) T = numpy.empty(Tmax-Tmin,numpy.int32) Y = numpy.empty(Tmax-Tmin) error = numpy.empty(Tmax-Tmin) dev = numpy.empty(Tmax-Tmin) like = numpy.empty(Tmax-Tmin) for t in xrange(Tmin,Tmax): fields = lines[t-Tmin+4].split() assert int(fields[0]) == t,"""first field should be %d found %s line='%s'"""%(t,fields[0],lines[t-Tmin+4]) T[t-Tmin] = int(fields[0]) Y[t-Tmin] = float(fields[1]) error[t-Tmin] = Y[t-Tmin]-float(fields[2]) dev[t-Tmin] = math.sqrt(float(fields[3])) like[t-Tmin] = float(fields[4]) return (T,Y,error,dev,like) def read_FU(name): """ Read a Forecast and Update from the file with name "name". The forecast is characterized by the mean "aM" and covariance "aSigma" and the update is characterized by mean "alphaM" and covariance "alphaSig" The file should have been written by Hview.py """ lines = open(name,'r').readlines() def skip_to(pattern): for i in xrange(len(lines)): if lines[i].startswith(pattern): return i raise RuntimeError,'%s not found'%pattern def read_matrix(pattern,Ni,Nj): import re start = skip_to(pattern) matrix = numpy.mat(numpy.empty((Ni,Nj))) for i in xrange(Ni): numbers = [float(s) for s in re.sub(r'[\]\[]*','', lines[start+1+i]).split()] for j in xrange(Nj): matrix[i,j] = numbers[j] return matrix # These try/except pairs let me read files from either the new or # old version of Hview try: alphaM = read_matrix('alphaM',3,1) aM = read_matrix('aM',3,1) except: alphaM = read_matrix('alphaM',1,3).T aM = read_matrix('aM',1,3).T try: alphaSig = read_matrix('alphaSig',3,3) except: alphaSig = INV(read_matrix('alphaSI',3,3)) aSigma = read_matrix('aSigma',3,3) return (aM,aSigma,alphaM,alphaSig) def ellipse(mu,cov): cov = numpy.mat(cov[::2,::2]) # drop component 1 mu = numpy.array(mu).reshape((3,))[::2] v,Q = EIG(cov) R = numpy.mat(numpy.diag(numpy.sqrt(v))) R = Q*R*INV(Q) # Square root of cov NC = 100.0 # Number of points on unit circle a = numpy.arange(0, 2*math.pi*(NC+1)/NC,2*math.pi/NC) C = numpy.mat(numpy.array([numpy.cos(a),numpy.sin(a)])) # A unit circle rv = (mu + (R*C).A.T).T assert rv.shape[0] == 2,'rv.shape=(%d,%d)'%rv.shape return rv def ellipse_FU(name): """ Project onto x_0,x_2. Calculate and return the following: ellipse_F Level set for aM, aSigma ellipse_U Level set for alphaM, alpha_sig x_center Average of min and max of first component of ellipses D_x Difference/2 between min and max of first component of ellipses y_center D_y """ aM,aSigma,alphaM,alphaSig = read_FU(name) ellipse_F = ellipse(aM,aSigma) ellipse_U = ellipse(alphaM+aM,alphaSig) maxes = numpy.maximum(ellipse_F.max(1),ellipse_U.max(1)) mins = numpy.minimum(ellipse_F.min(1),ellipse_U.min(1)) x_center,y_center = (maxes+mins)/2 D_x,D_y = (maxes-mins)/2 *1.2 return (ellipse_F,ellipse_U,x_center,D_x,y_center,D_y) T,Y,error,dev,like = read_TS(dataA) pylab.subplot(3,1,1) pylab.plot(T,Y,'g-') pylab.legend((r'$y(t)$',),loc='upper right') pylab.subplot(3,1,2) pylab.plot(T,dev,'r-',T,error,'b-') pylab.legend((r'$\sigma_\gamma(t)$',r'$y(t)-\mu_\gamma(t)$'),loc='upper right') pylab.subplot(3,1,3) pylab.plot(T,like,'m-') pylab.legend((r'$\log \left(P_\gamma(y(t))\right)$',),loc='lower right') pylab.xlabel(r'$t$') pylab.savefig(fig0) pylab.clf() # Clear figure window B_F,B_U,B_xc,B_Dx,B_yc,B_Dy = ellipse_FU(dataB) C_F,C_U,C_xc,C_Dx,C_yc,C_Dy = ellipse_FU(dataC) Dx = max(B_Dx,C_Dx) # X_width of both plots Dy = max(B_Dy,C_Dy) # Y_width of both plots pylab.subplot(1,2,1) pylab.plot(B_F[0],B_F[1],'c-',B_U[0],B_U[1],'r-') pylab.axis(xmin=B_xc-Dx,xmax=B_xc+Dx,ymin=B_yc-Dy,ymax=B_yc+Dy) pylab.legend(('Forecast','Update'),loc='upper left') pylab.subplot(1,2,2) pylab.plot(C_F[0],C_F[1],'r-',C_U[0],C_U[1],'b-') pylab.axis(xmin=C_xc-Dx,xmax=C_xc+Dx,ymin=C_yc-Dy,ymax=C_yc+Dy) pylab.legend(('Forecast','Update'),loc='upper center') pylab.savefig(fig1) # Local Variables: # mode: python # End: