In [1]:
from sklearn.datasets import fetch_olivetti_faces
from sklearn.preprocessing import StandardScaler
from sklearn.cross_validation import train_test_split
import OVFM.Model as md
import OVFM.FeatureMap as fm
import OVFM.Risk as rsk
import OVFM.LearningRate as lr
import OVFM.DataGeneration as dg
import OVFM.SGD as sgd
import numpy as np
import scipy as sp
import time as time
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import RidgeCV
In [93]:
def gaussian( X, Y, mean_X, mean_Y, scale = 1 ):
Xc = X - mean_X
Yc = Y - mean_Y
return np.pi ** 2 * np.exp( - scale / 2 * ( Xc ** 2 + Yc ** 2 ) ) / np.sqrt( scale )
def uniform_X( d, n ):
return 5 * ( 2 * np.random.rand( d, n ) - 1 )
def mesh2array( X, Y ):
arr = np.empty( ( X.size, 2 ) )
arr[ :, 0 ] = X.ravel( )
arr[ :, 1 ] = Y.ravel( )
return arr
def array2mesh( arr, side = None ):
if side == None:
side = int( np.sqrt( arr.shape[ 0 ] ) )
X = arr[ :, 0 ].reshape( ( side, side ) )
Y = arr[ :, 1 ].reshape( ( side, side ) )
return X, Y
def MSE( p, y ):
return np.mean( ( p - y ) ** 2 )
def MSE_DIR( p, y ):
return np.mean( ( ( p.T / np.linalg.norm( p, axis = 1 ) ).T - ( y.T / np.linalg.norm( y, axis = 1 ) ).T ) ** 2 )
def MSE_NORM( p, y ):
return np.mean( ( np.linalg.norm( p, axis = 1 ) - np.linalg.norm( y, axis = 1 ) ) ** 2 )
def RMSE( p, y ):
return np.sqrt( MSE( p, y ) )
def RMSE_DIR( p, y ):
return np.sqrt( MSE_DIR( p, y ) )
def RMSE_NORM( p, y ):
return np.sqrt( MSE_NORM( p, y ) )
In [142]:
# Generate data
d = 2
D = 1000
p = 2
d = 2
D = 1000
p = 2
x = np.arange( -1, 1, .01 )
y = np.arange( -1, 1, .01 )
X, Y = np.meshgrid( x, y )
loc = 40
space = 0.3
field = gaussian( X, Y, -space, 0, loc ) + gaussian( X, Y, space, 0, loc ) + gaussian( X, Y, 0, space, loc ) + gaussian( X, Y, 0, -space, loc )
V, U = np.gradient( field )
In [4]:
# Plot data
plt.streamplot( X, Y, U, V, color = np.sqrt( U ** 2 + V ** 2 ), linewidth = .5, cmap=plt.cm.autumn, density = 10, arrowstyle=u'->' )
plt.show( )
In [143]:
x = mesh2array( X, Y )
y = mesh2array( U, V )
In [144]:
A = 2 * np.random.rand( p, p ) - 1
A = np.dot( A, A.T )
A = A / np.linalg.norm( A, 2 )
yn = y + 0. * np.var( y ) * np.random.multivariate_normal( np.zeros( p ), A, y.shape[ 0 ] )
X_train, X_test, y_train, y_test = train_test_split( x, yn, train_size = .9 )
# X_train = x[ x[ :, 0 ] < 4.8, : ]
# y_train = y[ x[ :, 0 ] < 4.8, : ]
# X_test = x[ x[ :, 0 ] >= 4.8, : ]
# y_test = y[ x[ :, 0 ] >= 4.8, : ]
In [145]:
scaler_X = StandardScaler( )
X_train_s = scaler_X.fit_transform( X_train )
X_test_s = scaler_X.transform( X_test )
scaler_y = StandardScaler( )
y_train_s = scaler_y.fit_transform( y_train )
y_test_s = scaler_y.transform( y_test )
In [146]:
A = np.cov( y_train.T )
A = A / np.linalg.norm( A )
# A = np.eye( p )
sigma = np.median( sp.spatial.distance.squareform( sp.spatial.distance.pdist( X_train_s.T ) ) )
gamma = 1000000. / ( 2 * sigma ** 2 )
model = md.Model( fm.DecomposableFF( gamma, d, 500, A, eps_or_rank = 0 ) )
In [147]:
lbda1 = 0.00
lbda2 = 0.00
eps = 0.0
risk = rsk.ElasticNet( lbda1, lbda2, eps )
lb = lr.Constant( 0. ) # no bias
lc = lr.Constant( 1. )
opt = sgd.SGD( risk, 100000, lc, lb, 10, 0 )
model.reset( )
start = time.time( )
opt.fit( model, X_train_s.T, y_train.T )
stop = time.time( )
start_pred = time.clock( )
pred_train = model( X_train_s.T ).T
pred_test = model( X_test_s.T ).T
print lbda1, lbda2, RMSE( pred_train, y_train ), RMSE( pred_test, y_test ), stop - start
print lbda1, lbda2, RMSE_DIR( pred_train, y_train ), RMSE_DIR( pred_test, y_test ), stop - start
print lbda1, lbda2, RMSE_NORM( pred_train, y_train ), RMSE_NORM( pred_test, y_test ), stop - start
stop_pred = time.clock( )
print stop_pred - start_pred
In [148]:
sigma = np.median( sp.spatial.distance.squareform( sp.spatial.distance.pdist( X_train.T ) ) )
gamma = 1000000. / ( 2 * sigma ** 2 )
model = md.Model( fm.CurlFreeFF( gamma, d, 500 ) )
In [149]:
lbda1 = 0.00
lbda2 = 0.00
eps = 0.0
risk = rsk.ElasticNet( lbda1, lbda2, eps )
lb = lr.Constant( 0. ) # no bias
lc = lr.InvScaling( 1., lbda1, 1 )
opt = sgd.SGD( risk, 100000, lc, lb, 10, 0 )
model.reset( )
start = time.time( )
opt.fit( model, X_train.T, y_train.T )
stop = time.time( )
start_pred = time.clock( )
pred_train = model( X_train.T ).T
pred_test = model( X_test.T ).T
print lbda1, lbda2, RMSE( pred_train, y_train ), RMSE( pred_test, y_test ), stop - start
print lbda1, lbda2, RMSE_DIR( pred_train, y_train ), RMSE_DIR( pred_test, y_test ), stop - start
print lbda1, lbda2, RMSE_NORM( pred_train, y_train ), RMSE_NORM( pred_test, y_test ), stop - start
stop_pred = time.clock( )
print stop_pred - start_pred
In [58]:
UT, VT = array2mesh( model( mesh2array( X, Y ).T ).T )
In [138]:
plt.streamplot( X, Y, U, V, color = np.sqrt( UT ** 2 + VT ** 2 ), linewidth = 1, cmap=plt.cm.autumn, density = 10, arrowstyle=u'->' )
plt.show( )
In [64]:
plt.imshow( np.sqrt( ( UT - U ) ** 2 + ( VT - V ) ** 2 ) )
plt.show( )
In [204]:
d = 2
D = 1000
p = 2
x = np.arange( -1, 1, .01 )
y = np.arange( -1, 1, .01 )
X, Y = np.meshgrid( x, y )
loc = 40
space = 0.3
field = gaussian( X, Y, -space, 0, loc ) + gaussian( X, Y, space, 0, loc ) + gaussian( X, Y, 0, space, loc ) + gaussian( X, Y, 0, -space, loc )
V, U = np.gradient( field )
# plt.imshow( np.sqrt( ( U ) ** 2 + ( V ) ** 2 ) )
plt.imshow( field )
plt.show( )
In [107]:
plt.imshow( np.sqrt( ( U ) ** 2 + ( V ) ** 2 ) )
plt.show( )
In [105]:
gaussian( X, Y, 0, 0, loc ).max( )
Out[105]:
In [ ]: