In [4]:

    

import numpy as np

from bokeh.plotting import HBox, VBox, figure, show, output_file, GridPlot
from bokeh.models.mappers import LinearColorMapper
from bokeh.models import BasicTicker, Grid 

from sklearn import metrics
from sklearn import preprocessing
from sklearn.datasets import fetch_olivetti_faces
from sklearn.utils.validation import check_random_state
from sklearn.ensemble import ExtraTreesClassifier, RandomForestClassifier
from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge, RidgeCV
from sklearn.cross_validation import train_test_split
from sklearn.covariance import GraphLassoCV, ledoit_wolf
from sklearn.grid_search import GridSearchCV
from scipy.spatial import distance
 
import scipy
import sklearn
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 time
import sys


    

In [2]:

    

q_os = scipy.io.loadmat( 'GP_mag/data/q_os.mat' )
data = scipy.io.loadmat( 'GP_mag/data/scanOneDeskMagCorrected.mat' )


    

In [118]:

    

# from pylab import *
mag_field = data[ 'U' ][ 0, 0 ][ 3 ] # 0: acc, 1: gyro, 2: mag, 3: mag corrected
pos_field = data[ 'pos_mt_l' ]

N = mag_field.shape[ 0 ]
# pos_field = pos_field[ np.floor( N / 20 ):np.floor( 3.8 * N / 4 ), : ]
# mag_field = mag_field[ np.floor( N / 20 ):np.floor( 3.8 * N / 4 ), : ]
pos_field = pos_field[ np.floor( N / 10 ):np.floor( 3 * N / 4 ), : ]
mag_field = mag_field[ np.floor( N / 10 ):np.floor( 3 * N / 4 ), : ]

pos_field = preprocessing.StandardScaler( ).fit_transform( pos_field )
# mag_field = preprocessing.StandardScaler( ).fit_transform( mag_field )


    

In [119]:

    

output_file( "scanOneDeskMagCorrected.html", title = "scanOneDeskMagCorrected" )
p1 = figure( )
p1.line( pos_field[ :, 0 ], pos_field[ :, 1 ], line_color = "gray", line_width = 1, alpha = 0.2 )

speed = np.sqrt( mag_field[ :, 0 ] ** 2 + mag_field[ :, 1 ] ** 2 )
theta = np.arctan2( mag_field[ :, 1 ], mag_field[ :, 0 ] )
length = speed / 50
angle = theta
x1 = pos_field[ :, 0 ] + length * np.cos( angle )
y1 = pos_field[ :, 1 ] + length * np.sin( angle )

cm = np.array( [ "#C7E9B4", "#7FCDBB", "#41B6C4", "#1D91C0", "#225EA8", "#0C2C84" ] )
ix = ( ( length - length.min( ) ) / ( length.max( ) - length.min( ) ) * 5 ).astype( 'int' )
colors = cm[ ix ]

p1.segment( pos_field[ :, 0 ], pos_field[ :, 1 ], x1, y1, color=colors, line_width = 1 )

show( p1 )


    

In [81]:

    

samples = 10

phi_cf = fm.CurlFreeFF( 0.01, 3, 20 )
phi_df = fm.DivFreeFF( 0.01, 3, 20 )

phi_cf_pos_field = phi_cf.map( pos_field[::samples,:] )
phi_df_pos_field = phi_df.map( pos_field[::samples,:] )
phi_pos_field = np.asarray( np.bmat( [ [ phi_df_pos_field, np.zeros( phi_cf_pos_field.shape ) ], [ phi_df_pos_field, phi_cf_pos_field ] ] ) )


    

In [82]:

    

# phi_cf_pos_field_inv = np.linalg.pinv( phi_cf_pos_field )
# phi_df_pos_field_inv = np.linalg.pinv( phi_df_pos_field )
phi_pos_field_inv = np.linalg.pinv( phi_pos_field )


    

In [83]:

    

augmented_mag_field = np.concatenate( ( mag_field[::samples,:], np.zeros( mag_field[::samples,:].shape ) ), axis = 0 )
w = np.dot( phi_pos_field_inv, augmented_mag_field.ravel( ) )


    

In [84]:

    

augmented_mag_field_p = np.dot( phi_pos_field, w.T ).reshape( augmented_mag_field.shape )
print 'train MSE', ( ( augmented_mag_field - augmented_mag_field_p ) ** 2 ).mean( )
print 'train B-field MSE', 2 * ( ( augmented_mag_field[ :augmented_mag_field_p.shape [ 0 ] / 2, : ] - augmented_mag_field_p[ :augmented_mag_field_p.shape [ 0 ] / 2, : ] ) ** 2 ).mean( )
print 'train M-field MSE', 2 * ( ( augmented_mag_field[ augmented_mag_field_p.shape [ 0 ] / 2:, : ] - augmented_mag_field_p[ augmented_mag_field_p.shape [ 0 ] / 2:, : ] ) ** 2 ).mean( )


    

    




train MSE 0.0386053662855
train B-field MSE 0.143192074364
train M-field MSE 0.0112293907779

In [85]:

    

np.linalg.norm( w )


    

    
Out[85]:





1677849.3290207901

In [86]:

    

phi_cf_pos_field = phi_cf.map( pos_field[::samples/10,:] )
phi_df_pos_field = phi_df.map( pos_field[::samples/10,:] )
phi_pos_field = np.asarray( np.bmat( [ [ phi_df_pos_field, np.zeros( phi_cf_pos_field.shape ) ], [ phi_df_pos_field, phi_cf_pos_field ] ] ) )


    

In [87]:

    

augmented_mag_field = np.concatenate( ( mag_field[::samples/10,:], np.zeros( mag_field[::samples/10,:].shape ) ), axis = 0 )
augmented_mag_field_p = np.dot( phi_pos_field, w.T ).reshape( augmented_mag_field.shape )


    

In [88]:

    

print 'train MSE', ( ( augmented_mag_field - augmented_mag_field_p ) ** 2 ).mean( )
print 'train B-field MSE', 2 * ( ( augmented_mag_field[ :augmented_mag_field_p.shape [ 0 ] / 2, : ] - augmented_mag_field_p[ :augmented_mag_field_p.shape [ 0 ] / 2, : ] ) ** 2 ).mean( )
print 'train M-field MSE', 2 * ( ( augmented_mag_field[ augmented_mag_field_p.shape [ 0 ] / 2:, : ] - augmented_mag_field_p[ augmented_mag_field_p.shape [ 0 ] / 2:, : ] ) ** 2 ).mean( )


    

    




train MSE 0.0385874799476
train B-field MSE 0.143128086258
train M-field MSE 0.0112218335329

In [122]:

    

D = 1000
gamma = 1. / ( 2 * .15 ** 2 )
L = np.eye( 3 )
C = 0
eta0 = .05
samples = 10

df = fm.DecomposableFF( gamma, pos_field[::samples,:].shape[ 1 ], D, L )
dff = fm.DivFreeFF( gamma, pos_field[::samples,:].shape[ 1 ], D )
cff = fm.CurlFreeFF( gamma, pos_field[::samples,:].shape[ 1 ], D )
dcff = fm.MultipleFF( np.array( [ [ dff, None ], [ dff, cff ] ] ), np.array( [ .3, .3 ] ) )

modelD = md.Model( df )
modelDF = md.Model( dff )
modelCF = md.Model( cff )
modelCDF = md.Model( dcff )

print modelCDF.coefs.shape
risk = rsk.Ridge( C, 0 )
lc = lr.Constant( 1 * eta0 )
lb = lr.Constant( 0.0 * eta0 )
est = sgd.SGD( risk, 5.0, lc, lb, 10, 10000 )


    

    




(1, 6000)

In [123]:

    

modelCDF.reset( )
est.fit( modelCDF, pos_field[::samples,:], np.concatenate( ( mag_field[::samples,:], np.zeros( mag_field[::samples,:].shape ) ), axis = 1 ) )


    

    




0 0.176662008999 0.0 0.0 0.025033516202 0.0158957209727
10000 0.0980321570513 0.0 0.0 27.431461537 0.48447206738
20000 0.0861548525372 0.0 0.0 53.2974455713 0.794392945767
30000 0.0797886413953 0.0 0.0 77.4086116548 1.0315998154
40000 0.0753945607547 0.0 0.0 100.61588872 1.24753692479
50000 0.0723201787406 0.0 0.0 122.672752845 1.42233574628
60000 0.0698567693237 0.0 0.0 143.686981845 1.57280143429
70000 0.0677467401271 0.0 0.0 163.893582518 1.70940153973
80000 0.0657939159475 0.0 0.0 183.31815782 1.83818474482

In [125]:

    

print ( ( modelCDF( pos_field[::samples,:] )[:,:3] - mag_field[::samples,:] ) ** 2 ).mean( )
print ( modelCDF( pos_field[::samples,:] )[:,3:] ** 2 ).mean( )


    

    




0.0828378460191
0.0466684591538

In [126]:

    

print modelCDF( pos_field[::samples,:] )
print mag_field[::samples,:]


    

    




[[-0.03634714  0.34862037 -0.35537649  0.06083936 -0.01329143 -0.36864603]
 [-0.18303732  0.08892595 -0.49272757 -0.03759883 -0.28648934 -0.2592961 ]
 [-0.19826752 -0.11281998 -0.35859342  0.05654806 -0.12007922 -0.03179159]
 ..., 
 [-0.30332214  0.13927734 -0.31951174 -0.13966712 -0.01629391 -0.17133474]
 [-0.18549799  0.23797037 -0.17848419  0.04963258  0.11705744 -0.03321907]
 [-0.1093708   0.29615309 -0.11224422 -0.10249431 -0.12739036  0.21583442]]
[[-0.15284036 -0.19392526 -0.57789424]
 [-0.17368152 -0.25437289 -0.61890733]
 [-0.22255254 -0.31008217 -0.65476024]
 ..., 
 [-0.26870347  0.09228446 -0.21421789]
 [-0.28811503  0.21817608 -0.21839684]
 [-0.33955136  0.34169237 -0.34931831]]

In [281]:

    




    

    




1.00014410486

In [5]:

    

def quiver( title, X, Y, U, V, scale ):
    p = figure( title = title )

    speed = np.sqrt( U ** 2 + V ** 2 )
    theta = np.arctan2( V, U )
    length = speed * scale
    angle = theta
    X1 = X + length * np.cos( angle )
    Y1 = Y + length * np.sin( angle )

    cm = np.array( [ "#C7E9B4", "#7FCDBB", "#41B6C4", "#1D91C0", "#225EA8", "#0C2C84" ] )
    ix = ( ( length - length.min( ) ) / ( length.max( ) - length.min( ) ) * 5 ).ravel( ).astype( 'int' )
    colors = cm[ ix ]

    p.scatter( X.ravel( ), Y.ravel( ), size = 2 )
    p.segment( X.ravel( ), Y.ravel( ), X1.ravel( ), Y1.ravel( ), color=colors, line_width = 1 )
    return p


    

In [6]:

    

def magnetic_sphere( X, Y, a, M0 ):
    mu0 = 4 * np.pi * 1e-7

    rm = np.sqrt( X ** 2 + Y ** 2 )
    tm = np.arctan2( Y, X )
    
    MU = np.zeros( Y.shape )
    MV = np.full( X.shape, M0 )
    
    mU = 4. / 3. * np.pi * a ** 3 * MU * np.cos( tm )
    mV = 4. / 3. * np.pi * a ** 3 * MV * np.sin( tm )
    
    # Outside the sphere
    BU = mu0 / ( 4 * np.pi ) * ( -mU / rm ** 3 + 3 * ( mU * np.cos( tm ) + mV * np.sin( tm ) ) * np.cos( tm ) / rm ** 3 )
    BV = mu0 / ( 4 * np.pi ) * ( -mV / rm ** 3 + 3 * ( mU * np.cos( tm ) + mV * np.sin( tm ) ) * np.sin( tm ) / rm ** 3 )
    
    BU[ Y < 0 ] = -BU[ Y < 0 ]
    BV[ Y < 0 ] = -BV[ Y < 0 ]
    
    HU = BU / mu0
    HV = BV / mu0
    
    MU[ rm >= a ] = 0
    MV[ rm >= a ] = 0
    
    # Inside the sphere
    BU[ rm < a ] = ( 2. / 3. * mu0 * MU )[ rm < a ]
    BV[ rm < a ] = ( 2. / 3. * mu0 * MV )[ rm < a ]
    HU[ rm < a ] = ( -1. / 3. * MU )[ rm < a ]
    HV[ rm < a ] = ( -1. / 3. * MV )[ rm < a ]
    
    return MU, MV, BU, BV, HU, HV


    

In [7]:

    

x = np.linspace( -5, 5, 100 )
y = np.linspace( -5, 5, 100 )
Y, X = np.meshgrid( y, x )
mu0 = 4 * np.pi * 1e-7
a = 2.5
MU, MV, BU, BV, HU, HV = magnetic_sphere( X, Y, a, -1 )
BU = BU / mu0
BV = BV / mu0


    

In [283]:

    

output_file( "scanOneDeskMagCorrected.html", title = "scanOneDeskMagCorrected" )
p1 = quiver( 'M-field', X, Y, MU, MV, 0.5 )
p2 = quiver( 'B-field', X, Y, BU, BV, 0.5 )
p3 = quiver( 'H-field', X, Y, HU, HV, 0.5 )
show( HBox( p1, p2, p3 ) )


    

In [8]:

    

x = np.vstack( ( X.ravel( ), Y.ravel( ) ) ).T
y = np.vstack( ( BU.ravel( ), BV.ravel( ) ) ).T


    

In [9]:

    

x_out = x[ np.sqrt( x[ :, 0 ] ** 2 + x[ :, 1 ] ** 2 ) > a, : ]
y_out = y[ np.sqrt( x[ :, 0 ] ** 2 + x[ :, 1 ] ** 2 ) > a, : ]


    

In [10]:

    

D = 500
gamma = 1
L = np.eye( 2 )
C = 0
eta0 = 2.

df = fm.DecomposableFF( gamma, x_out.shape[ 1 ], D, L )
dff = fm.DivFreeFF( gamma, x_out.shape[ 1 ], D )
cff = fm.CurlFreeFF( gamma, x_out.shape[ 1 ], D )
dcff = fm.MultipleFF( np.array( [ [ dff, None ], [ dff, cff ] ] ), np.array( [ 1. / 4, 1. / 4 ] ) )

modelD = md.Model( df )
modelDF = md.Model( dff )
modelCF = md.Model( cff )
modelCDF = md.Model( dcff )
risk = rsk.Ridge( C, 0 )
lc = lr.Constant( 1. * eta0 )
lb = lr.Constant( 0.0 * eta0 )
est = sgd.SGD( risk, 1.0, lc, lb, 10, 10000 )


    

In [11]:

    

modelD.reset( )
est.fit( modelD, x_out, y_out )
yp = modelD( x )


    

    




0 0.00852512699506 0.0 0.0 0.0232332924599 0.00990592147648
10000 5.75906567029e-05 0.0 0.0 1.61441921826 0.175145886419
20000 3.37229657075e-05 0.0 0.0 1.86603166515 0.199917213458
30000 2.38060841364e-05 0.0 0.0 2.00756191575 0.212795066574
40000 1.81507487708e-05 0.0 0.0 2.14012520367 0.216711011483
50000 1.59129523683e-05 0.0 0.0 2.22001701321 0.217347616712
60000 1.42910099447e-05 0.0 0.0 2.32630676014 0.221000695504
70000 1.29172669627e-05 0.0 0.0 2.36767422344 0.219312436285
80000 1.22740085876e-05 0.0 0.0 2.40046524099 0.218182844653

In [12]:

    

BUDP = yp[ :, 0 ].reshape( ( 100, 100 ) )
BVDP = yp[ :, 1 ].reshape( ( 100, 100 ) )


    

In [13]:

    

modelDF.reset( )
est.fit( modelDF, x_out, y_out )
yp = modelDF( x )


    

    




0 0.0105316701829 0.0 0.0 0.0443663707549 0.0332185742738
10000 0.00163254215886 0.0 0.0 9.15471529013 0.701009524704
20000 0.0013900117394 0.0 0.0 11.4468710013 0.805009554671
30000 0.00131605833957 0.0 0.0 11.8659899152 0.821786797495
40000 0.00137055293093 0.0 0.0 12.2185839127 0.835983245036
50000 0.00138711348184 0.0 0.0 12.2671101692 0.841011529093
60000 0.00153573963111 0.0 0.0 12.3002697135 0.839549179257
70000 0.00140021711791 0.0 0.0 12.4301864715 0.84646160741
80000 0.00132919594567 0.0 0.0 12.3720654394 0.843749762569

In [14]:

    

BUDFP = yp[ :, 0 ].reshape( ( 100, 100 ) )
BVDFP = yp[ :, 1 ].reshape( ( 100, 100 ) )


    

In [15]:

    

modelCF.reset( )
est.fit( modelCF, x_out, y_out )
yp = modelCF( x )


    

    




0 0.0091882851128 0.0 0.0 0.0432180175541 0.0330041224418
10000 0.00022497600981 0.0 0.0 2.30454623023 0.212093825097
20000 0.000188824495746 0.0 0.0 2.69765311436 0.231074168392
30000 0.000181563701413 0.0 0.0 2.86563179412 0.234889795862
40000 0.000175769405191 0.0 0.0 2.98475582277 0.238949620886
50000 0.000175551530178 0.0 0.0 3.02516422953 0.23916681678
60000 0.000196617780487 0.0 0.0 3.13950500785 0.242057258098
70000 0.000172678480073 0.0 0.0 3.12568326072 0.236894276574
80000 0.000172183100975 0.0 0.0 3.18167750311 0.238517803421

In [16]:

    

BUCFP = yp[ :, 0 ].reshape( ( 100, 100 ) )
BVCFP = yp[ :, 1 ].reshape( ( 100, 100 ) )


    

In [17]:

    

modelCDF.reset( )
est.fit( modelCDF, x_out, np.concatenate( ( y_out, np.zeros( y_out.shape ) ), axis = 1 ) )
yp = modelCDF( x )


    

    




0 0.00510068541853 0.0 0.0 0.0443663707549 0.0332185742738
10000 0.00193362465396 0.0 0.0 42.6041933196 1.32339760885
20000 0.00110988213899 0.0 0.0 96.5751711431 2.0679744313
30000 0.000841758414705 0.0 0.0 137.100941917 2.50052585838
40000 0.000738200151015 0.0 0.0 165.844740858 2.78107149041
50000 0.000690800896041 0.0 0.0 185.134775396 2.95790317172
60000 0.000692801567092 0.0 0.0 197.990328946 3.0623369783
70000 0.000658805340916 0.0 0.0 208.553127585 3.15003757081
80000 0.00064603814209 0.0 0.0 214.070230119 3.19377347318

In [18]:

    

BUCDFP = yp[ :, 0 ].reshape( ( 100, 100 ) )
BVCDFP = yp[ :, 1 ].reshape( ( 100, 100 ) )

print yp[ :, 2: ]


    

    




[[-0.00493778  0.02895237]
 [-0.00014107  0.02593189]
 [ 0.00191307  0.02309701]
 ..., 
 [-0.00169648  0.02239162]
 [-0.0039959   0.02610208]
 [-0.00872924  0.03004249]]

In [19]:

    

XP = x[ :, 0 ].reshape( ( 100, 100 ) )
YP = x[ :, 1 ].reshape( ( 100, 100 ) )

out = XP ** 2 + YP ** 2 != a ** 2


    

In [20]:

    

output_file( "predicted_fields.html", title = "predicted_fields" )
p1 = quiver( 'B-field', XP[ out ], YP[ out ], BU[ out ], BV[ out ], 0.5 )
p1.line( np.linspace( -5, 5, 1000 ),  np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p1.line( np.linspace( -5, 5, 1000 ), -np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p2 = quiver( 'pred B-field decom', XP[ out ], YP[ out ], BUDP[ out ], BVDP[ out ], 0.5 )
p2.line( np.linspace( -5, 5, 1000 ),  np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p2.line( np.linspace( -5, 5, 1000 ), -np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p3 = quiver( 'pred B-field divf', XP[ out ], YP[ out ], BUDFP[ out ], BVDFP[ out ], 0.5 )
p3.line( np.linspace( -5, 5, 1000 ),  np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p3.line( np.linspace( -5, 5, 1000 ), -np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p4 = quiver( 'pred B-field curlf', XP[ out ], YP[ out ], BUCFP[ out ], BVCFP[ out ], 0.5 )
p4.line( np.linspace( -5, 5, 1000 ),  np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p4.line( np.linspace( -5, 5, 1000 ), -np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p5 = quiver( 'pred B-field divf-curlf', XP[ out ], YP[ out ], BUCDFP[ out ], BVCDFP[ out ], 0.5 )
p5.line( np.linspace( -5, 5, 1000 ),  np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
p5.line( np.linspace( -5, 5, 1000 ), -np.sqrt( a ** 2 - np.linspace( -5, 5, 1000 ) ** 2 ), color= "red" )
show( HBox( p1, p2, p3, p4, p5 ) )


    

In [32]:

    

A = np.random.rand( 5 )
print A
A = np.outer( A, A )
print A
print scipy.linalg.svd( A )


    

    




[ 0.50684537  0.56367471  0.30677107  0.54779019  0.87240125]
[[ 0.25689223  0.28569592  0.1554855   0.27764492  0.44217253]
 [ 0.28569592  0.31772918  0.1729191   0.30877548  0.49175053]
 [ 0.1554855   0.1729191   0.09410849  0.16804618  0.26762747]
 [ 0.27764492  0.30877548  0.16804618  0.30007409  0.47789285]
 [ 0.44217253  0.49175053  0.26762747  0.47789285  0.76108395]]
(array([[-0.38535988,  0.92031922,  0.03755192,  0.04940049, -0.02568531],
       [-0.42856783, -0.15097149, -0.25918918, -0.63126823, -0.57258939],
       [-0.23324128, -0.04962432,  0.11043986, -0.57673655,  0.77350754],
       [-0.41649066, -0.22099982,  0.85644178,  0.13222848, -0.16345554],
       [-0.66329587, -0.28092106, -0.43095398,  0.49894948,  0.21552267]]), array([  1.72988794e+00,   1.65981939e-16,   3.89161868e-17,
         3.13827449e-17,   4.51279288e-18]), array([[-0.38535988, -0.42856783, -0.23324128, -0.41649066, -0.66329587],
       [-0.57656679,  0.77546778, -0.1586853 , -0.20190816,  0.01650845],
       [-0.31024121, -0.21095611, -0.67166527,  0.61709506,  0.16524942],
       [-0.65008892, -0.33823842,  0.59197662,  0.11774913,  0.31413055],
       [ 0.01425795, -0.23679448, -0.34474184, -0.62537785,  0.65862029]]))

In [ ]: