<s> """<STR_LIT>""" <EOL> import numpy as np <EOL> import scipy . spatial . distance as ssd <EOL> import scipy . optimize as opt <EOL> from lfd . util import colorize <EOL> VERBOSE = False <EOL> ENABLE_SLOW_TESTS = False <EOL> def nan2zero ( x ) : <EOL> np . putmask ( x , np . isnan ( x ) , <NUM_LIT:0> ) <EOL> return x <EOL> def tps_apply_kernel ( distmat , dim ) : <EOL> """<STR_LIT>""" <EOL> if dim == <NUM_LIT:2> : <EOL> return <NUM_LIT:4> * distmat ** <NUM_LIT:2> * np . log ( distmat + <NUM_LIT> ) <EOL> elif dim == <NUM_LIT:3> : <EOL> return - distmat <EOL> else : <EOL> raise NotImplementedError <EOL> def tps_kernel_matrix ( x_na ) : <EOL> dim = x_na . shape [ <NUM_LIT:1> ] <EOL> distmat = ssd . squareform ( ssd . pdist ( x_na ) ) <EOL> return tps_apply_kernel ( distmat , dim ) <EOL> def tps_kernel_matrix2 ( x_na , y_ma ) : <EOL> dim = x_na . shape [ <NUM_LIT:1> ] <EOL> distmat = ssd . cdist ( x_na , y_ma ) <EOL> return tps_apply_kernel ( distmat , dim ) <EOL> def tps_eval ( x_ma , lin_ag , trans_g , w_ng , x_na ) : <EOL> K_mn = tps_kernel_matrix2 ( x_ma , x_na ) <EOL> return np . dot ( K_mn , w_ng ) + np . dot ( x_ma , lin_ag ) + trans_g [ None , : ] <EOL> def tps_grad ( x_ma , lin_ag , _trans_g , w_ng , x_na ) : <EOL> _N , D = x_na . shape <EOL> M = x_ma . shape [ <NUM_LIT:0> ] <EOL> assert x_ma . shape [ <NUM_LIT:1> ] == <NUM_LIT:3> <EOL> dist_mn = ssd . cdist ( x_ma , x_na , '<STR_LIT>' ) <EOL> grad_mga = np . empty ( ( M , D , D ) ) <EOL> lin_ga = lin_ag . T <EOL> for a in xrange ( D ) : <EOL> diffa_mn = x_ma [ : , a ] [ : , None ] - x_na [ : , a ] [ None , : ] <EOL> grad_mga [ : , : , a ] = lin_ga [ None , : , a ] - np . dot ( nan2zero ( diffa_mn / dist_mn ) , w_ng ) <EOL> return grad_mga <EOL> def tps_nr_grad ( x_ma , lin_ag , _trans_g , w_ng , x_na , return_tuple = False ) : <EOL> """<STR_LIT>""" <EOL> N , D = x_na . shape <EOL> M = x_ma . shape [ <NUM_LIT:0> ] <EOL> assert x_ma . shape [ <NUM_LIT:1> ] == <NUM_LIT:3> <EOL> dists_mn = ssd . cdist ( x_ma , x_na , '<STR_LIT>' ) <EOL> diffs_mna = x_ma [ : , None , : ] - x_na [ None , : , : ] <EOL> grad_mga = np . empty ( ( M , D , D ) ) <EOL> lin_ga = lin_ag . T <EOL> for a in xrange ( D ) : <EOL> grad_mga [ : , : , a ] = lin_ga [ None , : , a ] - np . dot ( nan2zero ( diffs_mna [ : , : , a ] / dists_mn ) , w_ng ) <EOL> Jw_mngab = - ( nan2zero ( diffs_mna [ : , : , None , : , None ] / dists_mn [ : , : , None , None , None ] ) ) * grad_mga [ : , None , : , None , : ] <EOL> Jw_mngab = Jw_mngab + Jw_mngab . transpose ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:2> , <NUM_LIT:4> , <NUM_LIT:3> ) <EOL> Jw_mabng = Jw_mngab . transpose ( <NUM_LIT:0> , <NUM_LIT:3> , <NUM_LIT:4> , <NUM_LIT:1> , <NUM_LIT:2> ) <EOL> Jw = Jw_mabng . reshape ( M * D ** <NUM_LIT:2> , N * D ) <EOL> Jl_mcgab = np . eye ( D ) [ None , : , None , : , None ] * grad_mga [ : , None , : , None , : ] <EOL> Jl_mcgab = Jl_mcgab + Jl_mcgab . transpose ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:2> , <NUM_LIT:4> , <NUM_LIT:3> ) <EOL> Jl_mabcg = Jl_mcgab . transpose ( <NUM_LIT:0> , <NUM_LIT:3> , <NUM_LIT:4> , <NUM_LIT:1> , <NUM_LIT:2> ) <EOL> Jl = Jl_mabcg . reshape ( M * D ** <NUM_LIT:2> , D * D ) <EOL> Jt = np . zeros ( ( M * D ** <NUM_LIT:2> , D ) ) <EOL> if return_tuple : <EOL> return Jl , Jt , Jw <EOL> else : <EOL> J = np . c_ [ Jl , Jt , Jw ] <EOL> return J <EOL> def tps_nr_err ( x_ma , lin_ag , trans_g , w_ng , x_na ) : <EOL> """<STR_LIT>""" <EOL> M , D = x_ma . shape <EOL> grad_mga = tps_grad ( x_ma , lin_ag , trans_g , w_ng , x_na ) <EOL> err_mab = np . empty ( ( M , D , D ) ) <EOL> for m in xrange ( M ) : <EOL> err_mab [ m ] = np . dot ( grad_mga [ m ] . T , grad_mga [ m ] ) - np . eye ( D ) <EOL> return err_mab . flatten ( ) <EOL> def tps_cost ( lin_ag , trans_g , w_ng , x_na , y_ng , bend_coef , K_nn = None , return_tuple = False , wt_n = None ) : <EOL> """<STR_LIT>""" <EOL> D = lin_ag . shape [ <NUM_LIT:0> ] <EOL> if K_nn is None : K_nn = tps_kernel_matrix ( x_na ) <EOL> if wt_n is None : wt_n = np . ones ( len ( x_na ) ) <EOL> ypred_ng = np . dot ( K_nn , w_ng ) + np . dot ( x_na , lin_ag ) + trans_g [ None , : ] <EOL> res_cost = ( wt_n [ : , None ] * ( ypred_ng - y_ng ) ** <NUM_LIT:2> ) . sum ( ) <EOL> bend_cost = bend_coef * sum ( np . dot ( w_ng [ : , g ] , np . dot ( K_nn , w_ng [ : , g ] ) ) for g in xrange ( D ) ) <EOL> if return_tuple : <EOL> return res_cost , bend_cost , res_cost + bend_cost <EOL> else : <EOL> return res_cost + bend_cost <EOL> def tps_nr_cost_eval ( lin_ag , trans_g , w_ng , x_na , y_ng , xnr_ma , bend_coef , nr_coef , K_nn = None , return_tuple = False ) : <EOL> D = lin_ag . shape [ <NUM_LIT:0> ] <EOL> if K_nn is None : K_nn = tps_kernel_matrix ( x_na ) <EOL> ypred_ng = np . dot ( K_nn , w_ng ) + np . dot ( x_na , lin_ag ) + trans_g [ None , : ] <EOL> res_cost = ( ( ypred_ng - y_ng ) ** <NUM_LIT:2> ) . sum ( ) <EOL> bend_cost = bend_coef * sum ( np . dot ( w_ng [ : , g ] , np . dot ( K_nn , w_ng [ : , g ] ) ) for g in xrange ( D ) ) <EOL> nr_cost = nr_coef * ( tps_nr_err ( xnr_ma , lin_ag , trans_g , w_ng , x_na ) ** <NUM_LIT:2> ) . sum ( ) <EOL> if return_tuple : <EOL> return res_cost , bend_cost , nr_cost , res_cost + bend_cost + nr_cost <EOL> else : <EOL> return res_cost + bend_cost + nr_cost <EOL> def tps_nr_cost_eval_general ( lin_ag , trans_g , w_eg , x_ea , y_ng , nr_ma , bend_coef , nr_coef , K_ee = None , return_tuple = True ) : <EOL> E , D = x_ea . shape <EOL> N = y_ng . shape [ <NUM_LIT:0> ] <EOL> M = nr_ma . shape [ <NUM_LIT:0> ] <EOL> assert E == N + <NUM_LIT:4> * M <EOL> K_ee = K_ee or tps_kernel_matrix ( x_ea ) <EOL> K_ne = K_ee [ : N ] <EOL> x_na = x_ea [ : N ] <EOL> ypred_ng = np . dot ( K_ne , w_eg ) + np . dot ( x_na , lin_ag ) + trans_g [ None , : ] <EOL> res_cost = ( ( ypred_ng - y_ng ) ** <NUM_LIT:2> ) . sum ( ) <EOL> bend_cost = bend_coef * sum ( np . dot ( w_eg [ : , g ] , np . dot ( - K_ee , w_eg [ : , g ] ) ) for g in xrange ( D ) ) <EOL> nr_cost = nr_coef * ( tps_nr_err ( nr_ma , lin_ag , trans_g , w_eg , x_ea ) ** <NUM_LIT:2> ) . sum ( ) <EOL> if return_tuple : <EOL> return res_cost , bend_cost , nr_cost , res_cost + bend_cost + nr_cost <EOL> else : <EOL> return res_cost + bend_cost + nr_cost <EOL> def tps_fit ( x_na , y_ng , bend_coef , rot_coef , wt_n = None , K_nn = None ) : <EOL> N , D = x_na . shape <EOL> K_nn = tps_kernel_matrix ( x_na ) if K_nn is None else K_nn <EOL> coef_ratio = bend_coef / rot_coef if rot_coef > <NUM_LIT:0> else <NUM_LIT:0> <EOL> A = np . zeros ( ( N + D + <NUM_LIT:1> , N + D + <NUM_LIT:1> ) ) <EOL> A [ : N , : N ] = K_nn <EOL> A . flat [ np . arange ( <NUM_LIT:0> , N ) * ( N + D + <NUM_LIT:2> ) ] += bend_coef / ( wt_n if wt_n is not None else <NUM_LIT:1> ) <EOL> A [ : N , N : N + D ] = x_na <EOL> A [ : N , N + D ] = <NUM_LIT:1> <EOL> A [ N : N + D , : N ] = x_na . T <EOL> A [ N + D , : N ] = <NUM_LIT:1> <EOL> A [ N : N + D , N : N + D ] = coef_ratio * np . eye ( D ) <EOL> B = np . empty ( ( N + D + <NUM_LIT:1> , D ) ) <EOL> B [ : N ] = y_ng <EOL> B [ N : N + D ] = coef_ratio * np . eye ( D ) <EOL> B [ N + D ] = <NUM_LIT:0> <EOL> X = np . linalg . solve ( A , B ) <EOL> w_ng = X [ : N , : ] <EOL> lin_ag = X [ N : N + D , : ] <EOL> trans_g = X [ N + D , : ] <EOL> return lin_ag , trans_g , w_ng <EOL> def solve_eqp1 ( H , f , A ) : <EOL> """<STR_LIT>""" <EOL> n_vars = H . shape [ <NUM_LIT:0> ] <EOL> assert H . shape [ <NUM_LIT:1> ] == n_vars <EOL> assert f . shape [ <NUM_LIT:0> ] == n_vars <EOL> assert A . shape [ <NUM_LIT:1> ] == n_vars <EOL> n_cnts = A . shape [ <NUM_LIT:0> ] <EOL> _u , _s , _vh = np . linalg . svd ( A . T ) <EOL> N = _u [ : , n_cnts : ] <EOL> L = N . T . dot ( H . dot ( N ) ) <EOL> R = - N . T . dot ( f ) <EOL> z = np . linalg . solve ( L , R ) <EOL> x = N . dot ( z ) <EOL> return x <EOL> def tps_fit3 ( x_na , y_ng , bend_coef , rot_coef , wt_n ) : <EOL> if wt_n is None : wt_n = np . ones ( len ( x_na ) ) <EOL> n , d = x_na . shape <EOL> K_nn = tps_kernel_matrix ( x_na ) <EOL> Q = np . c_ [ np . ones ( ( n , <NUM_LIT:1> ) ) , x_na , K_nn ] <EOL> rot_coefs = np . ones ( d ) * rot_coef if np . isscalar ( rot_coef ) else rot_coef <EOL> A = np . r_ [ np . zeros ( ( d + <NUM_LIT:1> , d + <NUM_LIT:1> ) ) , np . c_ [ np . ones ( ( n , <NUM_LIT:1> ) ) , x_na ] ] . T <EOL> solve_dim_separately = not np . isscalar ( bend_coef ) or ( wt_n . ndim > <NUM_LIT:1> and wt_n . shape [ <NUM_LIT:1> ] > <NUM_LIT:1> ) <EOL> if not solve_dim_separately : <EOL> WQ = wt_n [ : , None ] * Q <EOL> QWQ = Q . T . dot ( WQ ) <EOL> H = QWQ <EOL> H [ d + <NUM_LIT:1> : , d + <NUM_LIT:1> : ] += bend_coef * K_nn <EOL> H [ <NUM_LIT:1> : d + <NUM_LIT:1> , <NUM_LIT:1> : d + <NUM_LIT:1> ] += np . diag ( rot_coefs ) <EOL> f = - WQ . T . dot ( y_ng ) <EOL> f [ <NUM_LIT:1> : d + <NUM_LIT:1> , <NUM_LIT:0> : d ] -= np . diag ( rot_coefs ) <EOL> Theta = solve_eqp1 ( H , f , A ) <EOL> else : <EOL> bend_coefs = np . ones ( d ) * bend_coef if np . isscalar ( bend_coef ) else bend_coef <EOL> if wt_n . ndim == <NUM_LIT:1> : <EOL> wt_n = wt_n [ : , None ] <EOL> if wt_n . shape [ <NUM_LIT:1> ] == <NUM_LIT:1> : <EOL> wt_n = np . tile ( wt_n , ( <NUM_LIT:1> , d ) ) <EOL> Theta = np . empty ( ( <NUM_LIT:1> + d + n , d ) ) <EOL> for i in range ( d ) : <EOL> WQ = wt_n [ : , i ] [ : , None ] * Q <EOL> QWQ = Q . T . dot ( WQ ) <EOL> H = QWQ <EOL> H [ d + <NUM_LIT:1> : , d + <NUM_LIT:1> : ] += bend_coefs [ i ] * K_nn <EOL> H [ <NUM_LIT:1> : d + <NUM_LIT:1> , <NUM_LIT:1> : d + <NUM_LIT:1> ] += np . diag ( rot_coefs ) <EOL> f = - WQ . T . dot ( y_ng [ : , i ] ) <EOL> f [ <NUM_LIT:1> + i ] -= rot_coefs [ i ] <EOL> Theta [ : , i ] = solve_eqp1 ( H , f , A ) <EOL> return Theta [ <NUM_LIT:1> : d + <NUM_LIT:1> ] , Theta [ <NUM_LIT:0> ] , Theta [ d + <NUM_LIT:1> : ] <EOL> def tps_fit2 ( x_na , y_ng , bend_coef , rot_coef , wt_n = None ) : <EOL> if wt_n is not None : raise NotImplementedError <EOL> N , D = x_na . shape <EOL> _u , _s , _vh = np . linalg . svd ( np . c_ [ x_na , np . ones ( ( N , <NUM_LIT:1> ) ) ] , full_matrices = True ) <EOL> N_nq = _u [ : , <NUM_LIT:4> : ] <EOL> K_nn = tps_kernel_matrix ( x_na ) <EOL> Q_nn = np . c_ [ x_na , np . ones ( ( N , <NUM_LIT:1> ) ) , K_nn . dot ( N_nq ) ] <EOL> QQ_nn = np . dot ( Q_nn . T , Q_nn ) <EOL> A = QQ_nn <EOL> A [ <NUM_LIT:4> : , <NUM_LIT:4> : ] += bend_coef * N_nq . T . dot ( K_nn ) . dot ( N_nq ) <EOL> B = Q_nn . T . dot ( y_ng ) <EOL> A [ : <NUM_LIT:3> , : <NUM_LIT:3> ] += rot_coef * np . eye ( <NUM_LIT:3> ) <EOL> B [ : <NUM_LIT:3> , : <NUM_LIT:3> ] += rot_coef * np . eye ( <NUM_LIT:3> ) <EOL> X = np . linalg . solve ( A , B ) <EOL> lin_ag = X [ : D , : ] <EOL> trans_g = X [ D , : ] <EOL> w_ng = N_nq . dot ( X [ D + <NUM_LIT:1> : , : ] ) <EOL> return lin_ag , trans_g , w_ng <EOL> def tps_nr_fit ( x_na , y_ng , bend_coef , nr_ma , nr_coef , method = "<STR_LIT>" ) : <EOL> N , D = x_na . shape <EOL> lin_ag , trans_g , w_ng = tps_fit2 ( x_na , y_ng , bend_coef , <NUM_LIT> ) <EOL> _u , _s , _vh = np . linalg . svd ( np . c_ [ x_na , np . ones ( ( N , <NUM_LIT:1> ) ) ] , full_matrices = True ) <EOL> N_nq = _u [ : , <NUM_LIT:4> : ] <EOL> K_nn = tps_kernel_matrix ( x_na ) <EOL> Q_nn = np . c_ [ x_na , np . ones ( ( N , <NUM_LIT:1> ) ) , K_nn . dot ( N_nq ) ] <EOL> QQ_nn = np . dot ( Q_nn . T , Q_nn ) <EOL> Bend_nn = np . zeros ( ( N , N ) ) <EOL> Bend_nn [ <NUM_LIT:4> : , <NUM_LIT:4> : ] = - N_nq . T . dot ( K_nn . dot ( N_nq ) ) <EOL> n_iter = <NUM_LIT> <EOL> for i in xrange ( n_iter ) : <EOL> X_ng = np . r_ [ lin_ag , trans_g [ None , : ] , N_nq . T . dot ( w_ng ) ] <EOL> res_cost , bend_cost , nr_cost , fval = tps_nr_cost_eval ( lin_ag , trans_g , w_ng , x_na , y_ng , nr_ma , bend_coef , nr_coef , return_tuple = True ) <EOL> if VERBOSE : print colorize ( "<STR_LIT>" % ( i , fval ) , '<STR_LIT>' ) , <EOL> if VERBOSE : print "<STR_LIT>" % ( res_cost , bend_cost , nr_cost ) <EOL> Jl_zcg , Jt_zg , Jw_zng = tps_nr_grad ( nr_ma , lin_ag , trans_g , w_ng , x_na , return_tuple = True ) <EOL> nrerr_z = tps_nr_err ( nr_ma , lin_ag , trans_g , w_ng , x_na ) <EOL> if method == "<STR_LIT>" : <EOL> fullstep_ng = np . empty ( ( N , D ) ) <EOL> for g in xrange ( D ) : <EOL> J_zn = np . c_ [ Jl_zcg [ : , g : : D ] , Jt_zg [ : , g : : D ] , Jw_zng [ : , g : : D ] . dot ( N_nq ) ] <EOL> JJ_nn = np . dot ( J_zn . T , J_zn ) <EOL> A = nr_coef * JJ_nn + QQ_nn + bend_coef * Bend_nn <EOL> X0 = X_ng [ : , g ] <EOL> B = nr_coef * np . dot ( J_zn . T , np . dot ( J_zn , X0 ) - nrerr_z ) + Q_nn . T . dot ( y_ng [ : , g ] ) <EOL> fullstep_ng [ : , g ] = np . linalg . solve ( A , B ) - X_ng [ : , g ] <EOL> elif method == "<STR_LIT>" : <EOL> grad_ng = np . empty ( ( N , D ) ) <EOL> for g in xrange ( D - <NUM_LIT:1> , - <NUM_LIT:1> , - <NUM_LIT:1> ) : <EOL> Jnr_zn = np . c_ [ Jl_zcg [ : , g : : D ] , Jt_zg [ : , g : : D ] , Jw_zng [ : , g : : D ] . dot ( N_nq ) ] <EOL> grad_ng [ : , g ] = <NUM_LIT:2> * nr_coef * nrerr_z . dot ( Jnr_zn ) + <NUM_LIT:2> * Q_nn . T . dot ( Q_nn . dot ( X_ng [ : , g ] ) - y_ng [ : , g ] ) + <NUM_LIT:2> * bend_coef * Bend_nn . dot ( X_ng [ : , g ] ) <EOL> fullstep_ng = - grad_ng <EOL> cost_improved = False <EOL> for stepsize in <NUM_LIT> ** np . arange ( <NUM_LIT:0> , - <NUM_LIT:10> , - <NUM_LIT:1> ) : <EOL> cand_X_ng = X_ng + fullstep_ng * stepsize <EOL> cand_lin_ag , cand_trans_g , cand_w_ng = cand_X_ng [ : D ] , cand_X_ng [ D ] , N_nq . dot ( cand_X_ng [ D + <NUM_LIT:1> : ] ) <EOL> fval_cand = tps_nr_cost_eval ( cand_lin_ag , cand_trans_g , cand_w_ng , x_na , y_ng , nr_ma , bend_coef , nr_coef ) <EOL> if VERBOSE : print "<STR_LIT>" % ( stepsize , fval_cand ) <EOL> if fval_cand < fval : <EOL> cost_improved = True <EOL> break <EOL> if not cost_improved : <EOL> if VERBOSE : print "<STR_LIT>" <EOL> break <EOL> lin_ag = cand_lin_ag <EOL> trans_g = cand_trans_g <EOL> w_ng = cand_w_ng <EOL> return lin_ag , trans_g , w_ng <EOL> def tps_nr_fit_enhanced ( x_na , y_ng , bend_coef , nr_ma , nr_coef ) : <EOL> N , D = x_na . shape <EOL> M = nr_ma . shape [ <NUM_LIT:0> ] <EOL> E = N + <NUM_LIT:4> * M <EOL> F = E - M <EOL> Q = N + <NUM_LIT:3> * M - <NUM_LIT:4> <EOL> s = <NUM_LIT> <EOL> u = <NUM_LIT:1> / ( <NUM_LIT:2> * np . sqrt ( <NUM_LIT:2> ) ) <EOL> tetra_pts = [ ] <EOL> for pt in nr_ma : <EOL> tetra_pts . append ( s * np . r_ [ - <NUM_LIT> , <NUM_LIT:0> , - u ] + pt ) <EOL> tetra_pts . append ( s * np . r_ [ + <NUM_LIT> , <NUM_LIT:0> , - u ] + pt ) <EOL> tetra_pts . append ( s * np . r_ [ <NUM_LIT:0> , - <NUM_LIT> , + u ] + pt ) <EOL> tetra_pts . append ( s * np . r_ [ <NUM_LIT:0> , + <NUM_LIT> , + u ] + pt ) <EOL> x_ea = np . r_ [ x_na , tetra_pts ] <EOL> badsub_ex = np . c_ [ x_ea , np . ones ( ( E , <NUM_LIT:1> ) ) , np . r_ [ np . zeros ( ( N , M ) ) , np . repeat ( np . eye ( M ) , <NUM_LIT:4> , axis = <NUM_LIT:0> ) ] ] <EOL> lin_ag , trans_g , w_ng = tps_fit2 ( x_na , y_ng , bend_coef , <NUM_LIT> ) <EOL> w_eg = np . r_ [ w_ng , np . zeros ( ( <NUM_LIT:4> * M , D ) ) ] <EOL> assert badsub_ex . shape [ <NUM_LIT:0> ] >= badsub_ex . shape [ <NUM_LIT:1> ] <EOL> _u , _s , _vh = np . linalg . svd ( badsub_ex , full_matrices = True ) <EOL> assert badsub_ex . shape [ <NUM_LIT:1> ] == _s . size <EOL> N_eq = _u [ : , badsub_ex . shape [ <NUM_LIT:1> ] : ] <EOL> assert N_eq . shape == ( E , Q ) <EOL> assert E == N + <NUM_LIT:4> * M <EOL> assert F == Q + <NUM_LIT:4> <EOL> K_ee = tps_kernel_matrix ( x_ea ) <EOL> K_ne = K_ee [ : N , : ] <EOL> Q_nf = np . c_ [ x_na , np . ones ( ( N , <NUM_LIT:1> ) ) , K_ne . dot ( N_eq ) ] <EOL> QQ_ff = np . dot ( Q_nf . T , Q_nf ) <EOL> Bend_ff = np . zeros ( ( F , F ) ) <EOL> Bend_ff [ <NUM_LIT:4> : , <NUM_LIT:4> : ] = - N_eq . T . dot ( K_ee . dot ( N_eq ) ) <EOL> assert Q_nf . shape == ( N , F ) <EOL> assert w_eg . shape == ( E , D ) <EOL> n_iter = <NUM_LIT> <EOL> for i in xrange ( n_iter ) : <EOL> X_fg = np . r_ [ lin_ag , <EOL> trans_g [ None , : ] , <EOL> N_eq . T . dot ( w_eg ) ] <EOL> res_cost , bend_cost , nr_cost , fval = tps_nr_cost_eval_general ( lin_ag , trans_g , w_eg , x_ea , y_ng , nr_ma , bend_coef , nr_coef , return_tuple = True ) <EOL> if VERBOSE : print colorize ( "<STR_LIT>" % ( i , fval ) , '<STR_LIT>' ) , <EOL> if VERBOSE : print "<STR_LIT>" % ( res_cost , bend_cost , nr_cost ) <EOL> Jl_zcg , Jt_zg , Jw_zeg = tps_nr_grad ( nr_ma , lin_ag , trans_g , w_eg , x_ea , return_tuple = True ) <EOL> nrerr_z = tps_nr_err ( nr_ma , lin_ag , trans_g , w_eg , x_ea ) <EOL> fullstep_fg = np . empty ( ( F , D ) ) <EOL> for g in xrange ( D ) : <EOL> J_zf = np . c_ [ Jl_zcg [ : , g : : D ] , Jt_zg [ : , g : : D ] , Jw_zeg [ : , g : : D ] . dot ( N_eq ) ] <EOL> JJ_ff = np . dot ( J_zf . T , J_zf ) <EOL> A_ff = nr_coef * JJ_ff + QQ_ff + bend_coef * Bend_ff <EOL> X0 = X_fg [ : , g ] <EOL> B_f = nr_coef * np . dot ( J_zf . T , np . dot ( J_zf , X0 ) - nrerr_z ) + Q_nf . T . dot ( y_ng [ : , g ] ) <EOL> fullstep_fg [ : , g ] = np . linalg . solve ( A_ff , B_f ) - X_fg [ : , g ] <EOL> cost_improved = False <EOL> for stepsize in <NUM_LIT> ** np . arange ( <NUM_LIT:0> , - <NUM_LIT:10> , - <NUM_LIT:1> ) : <EOL> cand_X_fg = X_fg + fullstep_fg * stepsize <EOL> cand_lin_ag , cand_trans_g , cand_w_eg = cand_X_fg [ : D ] , cand_X_fg [ D ] , N_eq . dot ( cand_X_fg [ D + <NUM_LIT:1> : ] ) <EOL> fval_cand = tps_nr_cost_eval_general ( cand_lin_ag , cand_trans_g , cand_w_eg , x_ea , y_ng , nr_ma , bend_coef , nr_coef , return_tuple = False ) <EOL> if VERBOSE : print "<STR_LIT>" % ( stepsize , fval_cand ) <EOL> if fval_cand < fval : <EOL> cost_improved = True <EOL> break <EOL> if not cost_improved : <EOL> if VERBOSE : print "<STR_LIT>" <EOL> break <EOL> lin_ag = cand_lin_ag <EOL> trans_g = cand_trans_g <EOL> w_eg = cand_w_eg <EOL> return lin_ag , trans_g , w_eg , x_ea <EOL> def tps_fit_fixedrot ( x_na , y_ng , bend_coef , lin_ag , K_nn = None , wt_n = None ) : <EOL> """<STR_LIT>""" <EOL> if wt_n is not None : raise NotImplementedError <EOL> N , _D = x_na . shape <EOL> _u , _s , _vh = np . linalg . svd ( np . c_ [ x_na , np . ones ( ( N , <NUM_LIT:1> ) ) ] , full_matrices = True ) <EOL> N_nq = _u [ : , <NUM_LIT:4> : ] <EOL> K_nn = tps_kernel_matrix ( x_na ) <EOL> Q_nn = np . c_ [ np . ones ( ( N , <NUM_LIT:1> ) ) , K_nn . dot ( N_nq ) ] <EOL> QQ_nn = np . dot ( Q_nn . T , Q_nn ) <EOL> A = QQ_nn <EOL> A [ <NUM_LIT:1> : , <NUM_LIT:1> : ] += bend_coef * N_nq . T . dot ( K_nn ) . dot ( N_nq ) <EOL> B = Q_nn . T . dot ( y_ng - x_na . dot ( lin_ag ) ) <EOL> X = np . linalg . solve ( A , B ) <EOL> trans_g = X [ <NUM_LIT:0> , : ] <EOL> w_ng = N_nq . dot ( X [ <NUM_LIT:1> : , : ] ) <EOL> return trans_g , w_ng <EOL> def tps_fit_regrot ( x_na , y_ng , bend_coef , rfunc , wt_n = None , max_iter = <NUM_LIT:1> , inner_max_iter = <NUM_LIT:100> , rgrad = None , l_init = None ) : <EOL> """<STR_LIT>""" <EOL> if rgrad is not None : raise NotImplementedError <EOL> if wt_n is not None : raise NotImplementedError <EOL> N , _D = x_na . shape <EOL> K_nn = tps_kernel_matrix ( x_na ) <EOL> if l_init is None : <EOL> lin_ag , trans_g , w_ng = tps_fit3 ( x_na , y_ng , bend_coef , <NUM_LIT> , wt_n ) <EOL> else : <EOL> lin_ag = l_init <EOL> if True : print "<STR_LIT>" , lin_ag <EOL> trans_g , w_ng = tps_fit_fixedrot ( x_na , y_ng , bend_coef , lin_ag , K_nn , wt_n ) <EOL> Q_nn = np . eye ( N ) - np . outer ( np . ones ( N ) , np . ones ( N ) ) / N <EOL> for _ in xrange ( max_iter ) : <EOL> e_ng = y_ng - K_nn . dot ( w_ng ) <EOL> xQe_ag = x_na . T . dot ( Q_nn ) . dot ( e_ng ) <EOL> xQx_aa = x_na . T . dot ( Q_nn ) . dot ( x_na ) <EOL> eQe_g = ( ( e_ng - e_ng . mean ( axis = <NUM_LIT:0> ) ) ** <NUM_LIT:2> ) . sum ( axis = <NUM_LIT:0> ) <EOL> def f ( x ) : <EOL> b_ag = x . reshape ( <NUM_LIT:3> , <NUM_LIT:3> ) <EOL> out = sum ( [ b_ag [ : , i ] . T . dot ( xQx_aa ) . dot ( b_ag [ : , i ] ) - <NUM_LIT:2> * b_ag [ : , i ] . T . dot ( xQe_ag [ : , i ] ) + eQe_g [ i ] for i in xrange ( <NUM_LIT:3> ) ] ) + rfunc ( b_ag ) <EOL> return out <EOL> x0 = lin_ag . flatten ( ) <EOL> ( soln , fopt , _ , _ , _ , _ , allsolns ) = opt . fmin_powell ( f , x0 , maxiter = inner_max_iter , disp = bool ( VERBOSE ) , retall = True , full_output = True ) <EOL> if not np . isfinite ( soln ) . all ( ) : <EOL> print "<STR_LIT>" <EOL> soln = allsolns [ - <NUM_LIT:2> ] <EOL> if not np . isfinite ( fopt ) : <EOL> soln = lin_ag <EOL> print "<STR_LIT>" <EOL> lin_ag = soln . reshape ( <NUM_LIT:3> , <NUM_LIT:3> ) <EOL> assert np . isfinite ( lin_ag ) . all ( ) <EOL> trans_g , w_ng = tps_fit_fixedrot ( x_na , y_ng , bend_coef , lin_ag , K_nn , wt_n ) <EOL> if VERBOSE : print "<STR_LIT>" , lin_ag <EOL> return lin_ag , trans_g , w_ng <EOL> def tps_cost_regrot ( lin_ag , trans_g , w_ng , x_na , y_ng , bend_coef , rfunc , K_nn = None , wt_n = None ) : <EOL> """<STR_LIT>""" <EOL> if wt_n is not None : raise NotImplementedError <EOL> return tps_cost ( lin_ag , trans_g , w_ng , x_na , y_ng , bend_coef , K_nn ) + rfunc ( lin_ag ) </s>
