diff --git a/test/pade/padeOle.py b/test/pade/padeOle.py
new file mode 100644
--- /dev/null
+++ b/test/pade/padeOle.py
@@ -0,0 +1,294 @@
+#!/usr/bin/env python
+
+import professor2 as prof
+from numpy import linalg
+
+import numpy as np
+
+
+class Pade(object):
+ def __init__(self, X, Y, m, n):
+ self._Xraw=X
+ self._dim = len(X[0])
+ self._pmin=[self._Xraw[:,d].min() for d in range(self._dim)]
+ self._pmax=[self._Xraw[:,d].max() for d in range(self._dim)]
+ self._X=[]
+ for x in self._Xraw:
+ self._X.append([self.map_prange(x[d], self._pmin[d], self._pmax[d]) for d in range(self._dim)])
+
+ self._F=np.diag(Y)
+ self._m=m
+ self._n=n
+ self._K=m+n
+ self._M=prof.numCoeffs(self._dim, self._m)
+ self._N=prof.numCoeffs(self._dim, self._n)
+ self.mkMatrices()
+ self._struct_g = prof.mkStructure(self._dim, self._m)
+ self._struct_h = prof.mkStructure(self._dim, self._n)
+
+ def map_prange(self, x,a,b):
+ return (x-a)/float(b-a)
+
+
+ @property
+ def dim(self): return self._dim
+ @property
+ def M(self): return self._M
+ @property
+ def N(self): return self._N
+
+
+ def mkMatrices(self):
+ if self._dim==1:
+ self._VM=prof.mkProfessorMatrix([x for x in self._X], self._K)
+ else:
+ self._VM=prof.mkProfessorMatrix([x for x in self._X], self._K)
+ self._P=linalg.pinv(self._VM)
+ self._Z=self._P.dot(self._F).dot(self._VM)
+
+ self._Zhat = self._Z[ : , 0:(self._N)]
+ self._Zbar = self._Zhat[ 0:(self._M), : ]
+ if self._n==0:
+ self._Ztilde = self._Zhat[(self._M-1): , : ]
+ else:
+ self._Ztilde = self._Zhat[(self._M): , : ]
+
+ # from IPython import embed
+ # embed()
+
+ self._U,self._S,self._V = linalg.svd(self._Ztilde)
+
+ self._bcoeff = self._V[-1]
+ self._acoeff = self._Zbar.dot(self._bcoeff)
+
+
+ def val(self, Xraw):
+ X = [self.map_prange(Xraw[d], self._pmin[d], self._pmax[d]) for d in range(self._dim)]
+ lv_g = np.array(prof.mkLongVector(X, self._m, self._struct_g))
+ lv_h = np.array(prof.mkLongVector(X, self._n, self._struct_h))
+ g=self._acoeff.dot(lv_g)
+ h=self._bcoeff.dot(lv_h)
+ return g/h
+
+
+import sys
+m=int(sys.argv[2])
+n=int(sys.argv[3])
+D=np.loadtxt(sys.argv[1], dtype=np.float64)
+
+A = D[:,0:-1]
+Y = D[:,-1]
+
+import pylab
+
+# res=[]
+# for num, x in enumerate(A):
+ # res.append((PP.val(x)-Y[num])/Y[num])
+
+# pylab.hist(res)
+# pylab.savefig("pres.pdf")
+
+PP=Pade(A,Y,10,0)
+
+pylab.clf()
+pylab.plot(A,Y, "bx", label="Truth m=2 n=1")
+for mm in xrange(0,4):
+ for nn in xrange(0,3):
+ # for mm in xrange(1,5):
+ # :w
+ PP=Pade(A,Y,mm,nn)
+ PADE=np.array([PP.val(x) for x in A])
+ pylab.plot(A, PADE,label="P(m=%i n=%i)"%(mm, nn))
+
+
+pylab.xlabel("x")
+pylab.ylabel("f(x)")
+pylab.xscale("log")
+pylab.legend()
+pylab.savefig("cmppade.pdf")
+
+pylab.clf()
+pylab.plot(A,Y, "bx", label="Truth m=2 n=1")
+# PP=Pade(A,Y,2,1)
+# PADE=np.array([PP.val(x) for x in A])
+# pylab.plot(A, PADE,label="P(m=2 n=1)")
+
+X = A
+ipols = [prof.Ipol(X, Y, i, "order%i" % i) for i in [2,5,10]]
+Z = [np.array([i.value(x) for x in sorted(X)]) for i in ipols]
+for num, z in enumerate(Z):
+ pylab.plot(sorted(X), z, label="Poly order %i" % (ipols[num].order))
+
+
+PP=Pade(A,Y,10,0)
+PADE=np.array([PP.val(x) for x in A])
+pylab.plot(A, PADE, "k--", label="P(m=10 n=0)")
+
+pylab.xlabel("x")
+pylab.ylabel("f(x)")
+pylab.xscale("log")
+pylab.legend()
+pylab.savefig("cmpprof.pdf")
+
+sys.exit(1)
+
+# Legacy code
+
+K = m + n + 1
+
+
+
+dim = len(A[0])
+M=prof.numCoeffs(dim, m)
+N=prof.numCoeffs(dim, n)
+
+
+if len(A.shape)==1:
+ VM=prof.mkProfessorMatrix([[x] for x in A], K)
+else:
+ VM=prof.mkProfessorMatrix([x for x in A], K)
+
+# K=len(A)
+
+# def mkVanderMonde(params, K):
+ # import numpy as np
+ # V=np.zeros((len(params), K+1), dtype=np.float64)
+ # for i, p in enumerate(params):
+ # for j in xrange(K+1):
+ # # print i, j, V.shape
+ # V[i,j] = p**j
+ # return V
+
+# inputs=np.random.choice(len(A), len(A), replace=False)
+# print len(inputs), K
+# # VM = mkVanderMonde([A[x] for x in inputs], K)
+# VM = mkVanderMonde(A, K)
+
+F=np.diag(Y)
+# F=np.diag([Y[x] for x in inputs])
+# F=np.diag(Y[0:-1])
+# VM=VM[0:-1,:]
+
+from numpy import linalg
+# from scipy import linalg
+
+P=linalg.pinv(VM)
+
+Z=P.dot(F).dot(VM)
+
+# Zhat = Z[ : , 0:(n+1)]
+# Zbar = Zhat[ 0:m+1, : ]
+# Ztilde = Zhat[m+1: , : ]
+Zhat = Z[ : , 0:(N)]
+Zbar = Zhat[ 0:M, : ]
+Ztilde = Zhat[M: , : ]
+
+U,S,V = linalg.svd(Ztilde)
+# U,S,V = linalg.svd(Ztilde, lapack_driver='gesvd' )
+
+
+def f(x, A, B, m, n):
+ dim = len(A[0])
+ s_g = prof.mkStructure(dim, m)
+ l_g = prof.mkLongVector(x, m, s_g)
+
+
+ s_h = prof.mkStructure(dim, n)
+ l_h = prof.mkLongVector(x, n, s_h)
+
+
+ g=0
+ for numa, a in enumerate(A):
+ g+= a*x**numa
+ h=0
+ for numb, b in enumerate(B):
+ h+= b*x**numb
+ return g/h
+
+# for num, v in enumerate(V):
+ # print num
+
+ # bcoeff = v
+ # acoeff = Zbar.dot(bcoeff)
+ # print len(acoeff) + len(bcoeff), "coefficicents"
+
+ # import pylab
+ # pylab.clf()
+ # pylab.plot(A,Y, "bo", label="Da tru tru")
+ # pylab.plot(A,[f(x, acoeff,bcoeff) for x in A], "rx",label="Da pade %i"%num)
+
+ # pylab.xscale("log")
+
+ # pylab.legend()
+ # pylab.savefig("cmp_%i.pdf"%num)
+
+bcoeff = V[-1]
+acoeff = Zbar.dot(bcoeff)
+
+from IPython import embed
+embed()
+sys.exit(1)
+
+
+import numpy as np
+Y=np.array(Y)
+
+
+import pylab
+pylab.clf()
+pylab.plot(A,Y, "b.", label="True inputs")
+PADE=np.array([f(x, acoeff,bcoeff) for x in A])
+pylab.plot(A, PADE, "rx",label="Pade with %i coefficients"%ncoeff)
+
+X = np.array([[a] for a in A])
+ipols = [prof.Ipol(X, Y, i, "order%i" % i) for i in [2,5,10,15]]
+Z = [np.array([i.value(x) for x in sorted(X)]) for i in ipols]
+for num, z in enumerate(Z):
+ pylab.plot(sorted(X), z, label="Poly order %i" % (ipols[num].order))
+
+pylab.xlabel("x")
+pylab.ylabel("f(x)")
+pylab.xscale("log")
+# pylab.yscale("log")
+
+pylab.legend()
+pylab.savefig("cmpprof.pdf")
+
+pylab.clf()
+pylab.axhline(1)
+
+pylab.plot(A, PADE/Y, "rx",label="Pade with %i coefficients"%ncoeff)
+for num, z in enumerate(Z):
+ pylab.plot(sorted(X), z/Y, label="Ipol order %i" % (ipols[num].order))
+
+pylab.xlabel("x")
+pylab.ylabel("Ratio w.r.t. truth")
+
+pylab.xscale("log")
+pylab.ylim((-1,2))
+pylab.savefig("cmpratioprof.pdf")
+
+
+
+
+
+
+
+
+
+
+# from IPython import embed
+# embed()
+
+# Z = [[i.value(x) for x in sorted(X)] for i in ipols]
+
+# import pylab
+# pylab.plot(A, Y, "r.", label="Anchors")
+# for num, z in enumerate(Z):
+ # pylab.plot(sorted(X), z, label="Ipol order %i" % (ipols[num].order))
+# pylab.yscale(sys.argv[2])
+# pylab.xscale(sys.argv[3])
+# pylab.xlabel("M")
+# pylab.ylabel("f(M)")
+# pylab.legend(loc=1)
+# pylab.savefig(sys.argv[1].split(".")[0]+"_x_%s_y_%s.pdf"%(sys.argv[2], sys.argv[3]))