In [ ]:
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.utils import shuffle
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import RidgeCV
In [ ]:
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 [ ]:
# Generate data (one million points)
xs = np.arange( -1, 1, .002 )
ys = np.arange( -1, 1, .002 )
X, Y = np.meshgrid( xs, ys )
loc = 25
space = 0.5
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 [ ]:
fig = plt.figure()
plt.streamplot( X, Y, U, V, color = np.sqrt( U ** 2 + V ** 2 ), linewidth = .5, cmap=plt.cm.jet, density = 10, arrowstyle=u'->' )
fig.savefig( 'grad_field.eps', format='eps' )
In [ ]:
# Shuffle data, store it in x, y and compute sigma^2, gamma according to jaaco... heuristic
x = mesh2array( X, Y )
y = mesh2array( U, V )
xshuff, yshuff = shuffle( x, y )
sigma2 = np.percentile( sp.spatial.distance.pdist( xshuff[:1000,:], 'sqeuclidean' ), 1 )
gamma = 1. / ( 2. * sigma2 )
In [ ]:
# Add noise on the output and split randomly into train / validation / test set
d = 2
D = 1000
p = 2
A = 2 * np.random.rand( p, p ) - 1
A = np.dot( A, A.T )
A = A / np.linalg.norm( A, 2 )
yn = yshuff + 0. * np.std( y ) * np.random.multivariate_normal( np.zeros( p ), A, y.shape[ 0 ] )
X_train, X_test, y_train, y_test = train_test_split( xshuff, yn, train_size = .75 )
X_train, X_validation, y_train, y_validation = train_test_split( X_train, y_train, train_size = 0.666667 )
print X_train.shape
print X_validation.shape
print X_test.shape
In [ ]:
# gammas = [ 0.1 * gamma, gamma, 10 * gamma ]
gammas = np.logspace( np.log( gamma ) - 1, np.log( gamma ) + 1, 10, base = np.e )
lbdas = np.concatenate( ( np.zeros( 1 ), np.logspace( -10, 0, 10 ) ) )
eta0s = np.logspace( 3, -3, 10, base = np.e )
err = np.zeros( ( gammas.size, lbdas.size, eta0s.size, 3 ) )
for i, c_gamma in enumerate( gammas ):
A = np.cov( X_train.T )
A = A / np.linalg.norm( A, 2 )
models = [
md.Model( fm.DecomposableFF( c_gamma, d, D, A, eps_or_rank = 0 ) ),
md.Model( fm.DecomposableFF( c_gamma, d, D, np.eye( p ), eps_or_rank = 0 ) ),
md.Model( fm.CurlFreeFF( c_gamma, d, D ) ) ]
for j, c_lambda in enumerate( lbdas ):
for k, c_eta0 in enumerate( eta0s ):
for l, c_model in enumerate( models ):
risk = rsk.ElasticNet( c_lambda, 0, 0 )
lb = lr.Constant( 0. ) # no bias
lc = lr.InvScaling( c_eta0, c_lambda, 1 )
opt = sgd.SGD( risk, X_train.shape[ 0 ], lc, lb, 10, 0 )
c_model.reset( )
try:
opt.fit( c_model, X_train.T, y_train.T )
err[ i, j, k, l ] = RMSE( c_model( X_validation.T ).T, y_validation )
# fig = plt.figure()
# UP, VP = array2mesh( c_model( x.T ).T )
# if np.all( np.sqrt( UP ** 2 + VP ** 2 ) < np.inf ):
# plt.streamplot( X, Y, UP, VP, color = np.linalg.norm( np.dstack( [ UP, VP ] ), axis=-1 ), linewidth = .5, cmap=plt.cm.jet, density = 1, arrowstyle=u'->' )
# plt.xlim( xmax = 1, xmin = -1 )
# plt.ylim( ymax = 1, ymin = -1 )
# plt.colorbar( )
# fig.savefig( 'grad_field_' + str( i ) + '_' + str( j ) + '_' + str( k ) + '_' + str( l ) + '.eps', format='eps' )
except (RuntimeError):
err[ i, j, k, l ] = np.inf
# stop = time.time( )
In [ ]:
print UP
# plt.streamplot( X, Y, UP, VP, color = np.sqrt( UP ** 2 + VP ** 2 ), linewidth = .5, cmap=plt.cm.jet, density = 1, arrowstyle=u'->' )
In [ ]:
idx = 2
print np.unravel_index( np.nanargmin( err[ :, :, :, idx ].ravel( ) ), err[ :, :, :, idx ].shape )
print np.unravel_index( np.nanargmin( err[ :, :, :, : ].ravel( ) ), err[ :, :, :, : ].shape )
In [ ]:
print err[ 6, 1, 1, 0 ]
print err[ 6, 2, 1, 1 ]
print err[ 5, 1, 8, 2 ]
In [ ]:
print gammas
print lbdas
print eta0s
In [ ]:
opt.fit( c_model, X_train.T, y_train.T )
err[ i, j, k, l ] = RMSE( c_model( X_validation.T ).T, y_validation )
In [ ]:
print err
In [ ]:
print x
In [ ]: