Second pass at Ridge Opt

psdr/basis.py

@@ -320,8 +320,7 @@ def DDV(self, X):
 							eq = np.zeros(self.p+1)
 							eq[alpha[q]] = 1.
 							der2 = self.der(eq, 2)
-							#print "q", q, "der2", der2, eq
-							DDV[:,j,k,ell] *= V_coordinate[q][:,0:-2].dot(der2)
+							DDV[:,j,k,ell] *= V_coordinate[q][:,0:len(der2)].dot(der2)
 						elif q == k or q == ell:
 							DDV[:,j,k,ell] *= np.dot(V_coordinate[q][:,0:-1], self.Dmat[alpha[q],:])
 						else:

psdr/poly.py

@@ -44,15 +44,50 @@ def __init__(self, dimension, degree, coef):
 		self.basis = LegendreTensorBasis(dimension, degree) 
 		self.coef = coef
 
+	def V(self, X):	
+		return self.basis.V(X)
+
+	def DV(self, X):
+		return self.basis.DV(X)
+
+	def DDV(self, X):
+		return self.basis.DDV(X)
+
 	def eval(self, X):
-		V = self.basis.V(X)
-		return V.dot(self.coef) 
+		if len(X.shape) == 1:
+			return self.V(X.reshape(1,-1)).dot(self.coef).reshape(1)
+		else:
+			return self.V(X).dot(self.coef)
 
 	def grad(self, X):
-		return np.tensordot(self.basis.DV(X), self.coef, axes = (1,0))
+		if len(X.shape) == 1:
+			one_d = True
+			X = X.reshape(1,-1)	
+		else:
+			one_d = False	
+
+		DV = self.DV(X)
+		# Compute gradient on projected space
+		Df = np.tensordot(DV, self.coef, axes = (1,0))
+		# Inflate back to whole space
+		if one_d:
+			return Df.reshape(X.shape[1])
+		else:
+			return Df
 
 	def hessian(self, X):
-		return np.tensordot(self.basis.DDV(X), self.coef, axes = (1,0))
+		if len(X.shape) == 1:
+			one_d = True
+			X = X.reshape(1,-1)	
+		else:
+			one_d = False
+
+		DDV = self.DDV(X)
+		DDf = np.tensordot(DDV, self.coef, axes = (1,0))
+		if one_d:
+			return DDf.reshape(X.shape[1], X.shape[1])
+		else:
+			return DDf
 
 
 class PolynomialApproximation(PolynomialFunction):

psdr/polyridge.py

@@ -354,8 +354,15 @@ def fit(self, X, fX, U0 = None, **kwargs):
 	def _fit_affine(self, X, fX):
 		r""" Solves the affine 
 		"""
-		# TODO: There is often a scaling issue 
-		XX = np.hstack([X, np.ones((X.shape[0],1))])
+		# Normalize the domain 
+		lb = np.min(X, axis = 0)
+		ub = np.max(X, axis = 0)
+		dom = BoxDomain(lb, ub) 
+		XX = np.hstack([dom.normalize(X), np.ones((X.shape[0],1))])
+
+		# Normalize the output
+		fX = (fX - np.min(fX))/(np.max(fX) - np.min(fX))
+
 		if self.bound is None:
 			if self.norm == 1:
 				b = one_norm_fit(XX, fX)
@@ -370,6 +377,8 @@ def _fit_affine(self, X, fX):
 			b = bound_fit(-XX, -fX, norm = self.norm)	 	
 
 		U = b[0:-1].reshape(-1,1)
+		# Correct for transform 
+		U = dom._normalize_der().dot(U)
 		return U