<s> import copy <EOL> import warnings <EOL> from math import sqrt , exp , log , cosh , sinh <EOL> import numpy as np <EOL> from scipy import linalg <EOL> from numpy . testing import assert_array_almost_equal , assert_array_equal <EOL> from nose . tools import assert_raises , assert_equal , assert_true <EOL> from sklearn . utils import check_random_state <EOL> from sklearn . covariance import EmpiricalCovariance , LedoitWolf <EOL> from nilearn . _utils . extmath import is_spd <EOL> from nilearn . connectome . connectivity_matrices import ( <EOL> _check_square , _check_spd , _map_eigenvalues , _form_symmetric , <EOL> _geometric_mean , sym_to_vec , _prec_to_partial , ConnectivityMeasure ) <EOL> def grad_geometric_mean ( mats , init = None , max_iter = <NUM_LIT:10> , tol = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> mats = np . array ( mats ) <EOL> if init is None : <EOL> gmean = np . mean ( mats , axis = <NUM_LIT:0> ) <EOL> else : <EOL> gmean = init <EOL> norm_old = np . inf <EOL> step = <NUM_LIT:1.> <EOL> grad_norm = [ ] <EOL> for n in range ( max_iter ) : <EOL> vals_gmean , vecs_gmean = linalg . eigh ( gmean ) <EOL> gmean_inv_sqrt = _form_symmetric ( np . sqrt , <NUM_LIT:1.> / vals_gmean , vecs_gmean ) <EOL> whitened_mats = [ gmean_inv_sqrt . dot ( mat ) . dot ( gmean_inv_sqrt ) <EOL> for mat in mats ] <EOL> logs = [ _map_eigenvalues ( np . log , w_mat ) for w_mat in whitened_mats ] <EOL> logs_mean = np . mean ( logs , axis = <NUM_LIT:0> ) <EOL> norm = np . linalg . norm ( logs_mean ) <EOL> vals_log , vecs_log = linalg . eigh ( logs_mean ) <EOL> gmean_sqrt = _form_symmetric ( np . sqrt , vals_gmean , vecs_gmean ) <EOL> gmean = gmean_sqrt . dot ( <EOL> _form_symmetric ( np . exp , vals_log * step , vecs_log ) ) . dot ( gmean_sqrt ) <EOL> if norm < norm_old : <EOL> norm_old = norm <EOL> if norm > norm_old : <EOL> step = step / <NUM_LIT> <EOL> norm = norm_old <EOL> grad_norm . append ( norm / gmean . size ) <EOL> if tol is not None and norm / gmean . size < tol : <EOL> break <EOL> return grad_norm <EOL> def test_check_square ( ) : <EOL> non_square = np . ones ( ( <NUM_LIT:2> , <NUM_LIT:3> ) ) <EOL> assert_raises ( ValueError , _check_square , non_square ) <EOL> def test_check_spd ( ) : <EOL> non_sym = np . array ( [ [ <NUM_LIT:0> , <NUM_LIT:1> ] , [ <NUM_LIT:0> , <NUM_LIT:0> ] ] ) <EOL> assert_raises ( ValueError , _check_spd , non_sym ) <EOL> non_spd = np . ones ( ( <NUM_LIT:3> , <NUM_LIT:3> ) ) <EOL> assert_raises ( ValueError , _check_spd , non_spd ) <EOL> def test_map_eigenvalues ( ) : <EOL> sym = np . ones ( ( <NUM_LIT:2> , <NUM_LIT:2> ) ) <EOL> sym_exp = exp ( <NUM_LIT:1.> ) * np . array ( [ [ cosh ( <NUM_LIT:1.> ) , sinh ( <NUM_LIT:1.> ) ] , [ sinh ( <NUM_LIT:1.> ) , cosh ( <NUM_LIT:1.> ) ] ] ) <EOL> assert_array_almost_equal ( _map_eigenvalues ( np . exp , sym ) , sym_exp ) <EOL> spd_sqrt = np . array ( [ [ <NUM_LIT> , - <NUM_LIT:1.> , <NUM_LIT:0.> ] , [ - <NUM_LIT:1.> , <NUM_LIT> , - <NUM_LIT:1.> ] , [ <NUM_LIT:0.> , - <NUM_LIT:1.> , <NUM_LIT> ] ] ) <EOL> spd = spd_sqrt . dot ( spd_sqrt ) <EOL> assert_array_almost_equal ( _map_eigenvalues ( np . sqrt , spd ) , spd_sqrt ) <EOL> spd = np . array ( [ [ <NUM_LIT> , <NUM_LIT> ] , [ <NUM_LIT> , <NUM_LIT> ] ] ) <EOL> spd_log = np . array ( [ [ <NUM_LIT:0.> , log ( <NUM_LIT> ) ] , [ log ( <NUM_LIT> ) , <NUM_LIT:0.> ] ] ) <EOL> assert_array_almost_equal ( _map_eigenvalues ( np . log , spd ) , spd_log ) <EOL> def test_geometric_mean_couple ( ) : <EOL> n_features = <NUM_LIT:7> <EOL> spd1 = np . ones ( ( n_features , n_features ) ) <EOL> spd1 = spd1 . dot ( spd1 ) + n_features * np . eye ( n_features ) <EOL> spd2 = np . tril ( np . ones ( ( n_features , n_features ) ) ) <EOL> spd2 = spd2 . dot ( spd2 . T ) <EOL> vals_spd2 , vecs_spd2 = np . linalg . eigh ( spd2 ) <EOL> spd2_sqrt = _form_symmetric ( np . sqrt , vals_spd2 , vecs_spd2 ) <EOL> spd2_inv_sqrt = _form_symmetric ( np . sqrt , <NUM_LIT:1.> / vals_spd2 , vecs_spd2 ) <EOL> geo = spd2_sqrt . dot ( _map_eigenvalues ( np . sqrt , spd2_inv_sqrt . dot ( spd1 ) . dot ( <EOL> spd2_inv_sqrt ) ) ) . dot ( spd2_sqrt ) <EOL> assert_array_almost_equal ( _geometric_mean ( [ spd1 , spd2 ] ) , geo ) <EOL> def test_geometric_mean_diagonal ( ) : <EOL> n_matrices = <NUM_LIT:20> <EOL> n_features = <NUM_LIT:5> <EOL> diags = [ ] <EOL> for k in range ( n_matrices ) : <EOL> diag = np . eye ( n_features ) <EOL> diag [ k % n_features , k % n_features ] = <NUM_LIT> + k <EOL> diag [ ( n_features - <NUM_LIT:1> ) // ( k + <NUM_LIT:1> ) , ( n_features - <NUM_LIT:1> ) // ( k + <NUM_LIT:1> ) ] = ( k + <NUM_LIT:1> ) * <NUM_LIT> <EOL> diags . append ( diag ) <EOL> geo = np . prod ( np . array ( diags ) , axis = <NUM_LIT:0> ) ** ( <NUM_LIT:1> / float ( len ( diags ) ) ) <EOL> assert_array_almost_equal ( _geometric_mean ( diags ) , geo ) <EOL> def test_geometric_mean_geodesic ( ) : <EOL> n_matrices = <NUM_LIT:10> <EOL> n_features = <NUM_LIT:6> <EOL> sym = np . arange ( n_features ) / np . linalg . norm ( np . arange ( n_features ) ) <EOL> sym = sym * sym [ : , np . newaxis ] <EOL> times = np . arange ( n_matrices ) <EOL> non_singular = np . eye ( n_features ) <EOL> non_singular [ <NUM_LIT:1> : <NUM_LIT:3> , <NUM_LIT:1> : <NUM_LIT:3> ] = np . array ( [ [ - <NUM_LIT:1> , - <NUM_LIT> ] , [ - <NUM_LIT> , - <NUM_LIT:1> ] ] ) <EOL> spds = [ ] <EOL> for time in times : <EOL> spds . append ( non_singular . dot ( _map_eigenvalues ( np . exp , time * sym ) ) . dot ( <EOL> non_singular . T ) ) <EOL> gmean = non_singular . dot ( _map_eigenvalues ( np . exp , times . mean ( ) * sym ) ) . dot ( <EOL> non_singular . T ) <EOL> assert_array_almost_equal ( _geometric_mean ( spds ) , gmean ) <EOL> def random_diagonal ( p , v_min = <NUM_LIT:1.> , v_max = <NUM_LIT> , random_state = <NUM_LIT:0> ) : <EOL> """<STR_LIT>""" <EOL> random_state = check_random_state ( random_state ) <EOL> diag = random_state . rand ( p ) * ( v_max - v_min ) + v_min <EOL> diag [ diag == np . amax ( diag ) ] = v_max <EOL> diag [ diag == np . amin ( diag ) ] = v_min <EOL> return np . diag ( diag ) <EOL> def random_spd ( p , eig_min , cond , random_state = <NUM_LIT:0> ) : <EOL> """<STR_LIT>""" <EOL> random_state = check_random_state ( random_state ) <EOL> mat = random_state . randn ( p , p ) <EOL> unitary , _ = linalg . qr ( mat ) <EOL> diag = random_diagonal ( p , v_min = eig_min , v_max = cond * eig_min , <EOL> random_state = random_state ) <EOL> return unitary . dot ( diag ) . dot ( unitary . T ) <EOL> def random_non_singular ( p , sing_min = <NUM_LIT:1.> , sing_max = <NUM_LIT> , random_state = <NUM_LIT:0> ) : <EOL> """<STR_LIT>""" <EOL> random_state = check_random_state ( random_state ) <EOL> diag = random_diagonal ( p , v_min = sing_min , v_max = sing_max , <EOL> random_state = random_state ) <EOL> mat1 = random_state . randn ( p , p ) <EOL> mat2 = random_state . randn ( p , p ) <EOL> unitary1 , _ = linalg . qr ( mat1 ) <EOL> unitary2 , _ = linalg . qr ( mat2 ) <EOL> return unitary1 . dot ( diag ) . dot ( unitary2 . T ) <EOL> def test_geometric_mean_properties ( ) : <EOL> n_matrices = <NUM_LIT> <EOL> n_features = <NUM_LIT:15> <EOL> spds = [ ] <EOL> for k in range ( n_matrices ) : <EOL> spds . append ( random_spd ( n_features , eig_min = <NUM_LIT:1.> , cond = <NUM_LIT> , <EOL> random_state = <NUM_LIT:0> ) ) <EOL> input_spds = copy . copy ( spds ) <EOL> gmean = _geometric_mean ( spds ) <EOL> assert_true ( isinstance ( spds , list ) ) <EOL> for spd , input_spd in zip ( spds , input_spds ) : <EOL> assert_array_equal ( spd , input_spd ) <EOL> assert ( is_spd ( gmean , decimal = <NUM_LIT:7> ) ) <EOL> spds . reverse ( ) <EOL> spds . insert ( <NUM_LIT:0> , spds [ <NUM_LIT:1> ] ) <EOL> spds . pop ( <NUM_LIT:2> ) <EOL> assert_array_almost_equal ( _geometric_mean ( spds ) , gmean ) <EOL> non_singular = random_non_singular ( n_features , random_state = <NUM_LIT:0> ) <EOL> spds_cong = [ non_singular . dot ( spd ) . dot ( non_singular . T ) for spd in spds ] <EOL> assert_array_almost_equal ( _geometric_mean ( spds_cong ) , <EOL> non_singular . dot ( gmean ) . dot ( non_singular . T ) ) <EOL> spds_inv = [ linalg . inv ( spd ) for spd in spds ] <EOL> init = linalg . inv ( np . mean ( spds , axis = <NUM_LIT:0> ) ) <EOL> assert_array_almost_equal ( _geometric_mean ( spds_inv , init = init ) , <EOL> linalg . inv ( gmean ) ) <EOL> grad_norm = grad_geometric_mean ( spds , tol = <NUM_LIT> ) <EOL> difference = np . diff ( grad_norm ) <EOL> assert_true ( np . amax ( difference ) <= <NUM_LIT:0.> ) <EOL> max_iter = <NUM_LIT:1> <EOL> tol = <NUM_LIT> <EOL> with warnings . catch_warnings ( record = True ) as w : <EOL> warnings . simplefilter ( "<STR_LIT>" ) <EOL> gmean = _geometric_mean ( spds , max_iter = max_iter , tol = tol ) <EOL> assert_equal ( len ( w ) , <NUM_LIT:1> ) <EOL> grad_norm = grad_geometric_mean ( spds , max_iter = max_iter , tol = tol ) <EOL> assert_equal ( len ( grad_norm ) , max_iter ) <EOL> assert_true ( grad_norm [ - <NUM_LIT:1> ] > tol ) <EOL> for p in [ <NUM_LIT> , <NUM_LIT:1.> ] : <EOL> spds = [ ] <EOL> for k in range ( int ( p * n_matrices ) ) : <EOL> spds . append ( random_spd ( n_features , eig_min = <NUM_LIT> , cond = <NUM_LIT> , <EOL> random_state = <NUM_LIT:0> ) ) <EOL> for k in range ( int ( p * n_matrices ) , n_matrices ) : <EOL> spds . append ( random_spd ( n_features , eig_min = <NUM_LIT:1.> , cond = <NUM_LIT> , <EOL> random_state = <NUM_LIT:0> ) ) <EOL> if p < <NUM_LIT:1> : <EOL> max_iter = <NUM_LIT:30> <EOL> else : <EOL> max_iter = <NUM_LIT> <EOL> gmean = _geometric_mean ( spds , max_iter = max_iter , tol = <NUM_LIT> ) <EOL> def test_geometric_mean_errors ( ) : <EOL> n_features = <NUM_LIT:5> <EOL> mat1 = np . ones ( ( n_features , n_features + <NUM_LIT:1> ) ) <EOL> assert_raises ( ValueError , _geometric_mean , [ mat1 ] ) <EOL> mat1 = np . eye ( n_features ) <EOL> mat2 = np . ones ( ( n_features + <NUM_LIT:1> , n_features + <NUM_LIT:1> ) ) <EOL> assert_raises ( ValueError , _geometric_mean , [ mat1 , mat2 ] ) <EOL> assert_raises ( ValueError , _geometric_mean , [ mat2 ] ) <EOL> def test_sym_to_vec ( ) : <EOL> sym = np . ones ( ( <NUM_LIT:3> , <NUM_LIT:3> ) ) <EOL> vec = np . array ( [ <NUM_LIT:1.> , sqrt ( <NUM_LIT:2> ) , <NUM_LIT:1.> , sqrt ( <NUM_LIT:2> ) , sqrt ( <NUM_LIT:2> ) , <NUM_LIT:1.> ] ) <EOL> assert_array_almost_equal ( sym_to_vec ( sym ) , vec ) <EOL> def test_prec_to_partial ( ) : <EOL> prec = np . array ( [ [ <NUM_LIT> , - <NUM_LIT:1.> , <NUM_LIT:1.> ] , [ - <NUM_LIT:1.> , <NUM_LIT> , - <NUM_LIT:1.> ] , [ <NUM_LIT:1.> , - <NUM_LIT:1.> , <NUM_LIT:1.> ] ] ) <EOL> partial = np . array ( [ [ <NUM_LIT:1.> , <NUM_LIT> , - sqrt ( <NUM_LIT> ) / <NUM_LIT> ] , [ <NUM_LIT> , <NUM_LIT:1.> , sqrt ( <NUM_LIT> ) / <NUM_LIT> ] , <EOL> [ - sqrt ( <NUM_LIT> ) / <NUM_LIT> , sqrt ( <NUM_LIT> ) / <NUM_LIT> , <NUM_LIT:1.> ] ] ) <EOL> assert_array_almost_equal ( _prec_to_partial ( prec ) , partial ) <EOL> def test_connectivity_measure_errors ( ) : <EOL> conn_measure = ConnectivityMeasure ( ) <EOL> assert_raises ( ValueError , conn_measure . fit , <NUM_LIT:1.> ) <EOL> assert_raises ( ValueError , conn_measure . fit , [ np . ones ( ( <NUM_LIT:100> , <NUM_LIT> ) ) , <EOL> np . ones ( ( <NUM_LIT:10> , ) ) ] ) <EOL> assert_raises ( ValueError , conn_measure . fit , <EOL> [ np . ones ( ( <NUM_LIT:100> , <NUM_LIT> ) ) , np . ones ( ( <NUM_LIT:100> , <NUM_LIT> ) ) ] ) <EOL> def test_connectivity_measure_outputs ( ) : <EOL> n_subjects = <NUM_LIT:10> <EOL> n_features = <NUM_LIT> <EOL> n_samples = <NUM_LIT:200> <EOL> emp_covs = [ ] <EOL> ledoit_covs = [ ] <EOL> signals = [ ] <EOL> random_state = check_random_state ( <NUM_LIT:0> ) <EOL> ledoit_estimator = LedoitWolf ( ) <EOL> for k in range ( n_subjects ) : <EOL> signal = random_state . randn ( n_samples , n_features ) <EOL> signals . append ( signal ) <EOL> signal -= signal . mean ( axis = <NUM_LIT:0> ) <EOL> emp_covs . append ( ( signal . T ) . dot ( signal ) / n_samples ) <EOL> ledoit_covs . append ( ledoit_estimator . fit ( signal ) . covariance_ ) <EOL> kinds = [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" ] <EOL> for cov_estimator , covs in zip ( [ EmpiricalCovariance ( ) , LedoitWolf ( ) ] , <EOL> [ emp_covs , ledoit_covs ] ) : <EOL> input_covs = copy . copy ( covs ) <EOL> for kind in kinds : <EOL> conn_measure = ConnectivityMeasure ( kind = kind , <EOL> cov_estimator = cov_estimator ) <EOL> connectivities = conn_measure . fit_transform ( signals ) <EOL> assert_true ( isinstance ( connectivities , np . ndarray ) ) <EOL> assert_equal ( len ( connectivities ) , len ( covs ) ) <EOL> for k , cov_new in enumerate ( connectivities ) : <EOL> assert_array_equal ( input_covs [ k ] , covs [ k ] ) <EOL> assert ( is_spd ( covs [ k ] , decimal = <NUM_LIT:7> ) ) <EOL> if kind == "<STR_LIT>" : <EOL> assert_array_almost_equal ( cov_new , cov_new . T ) <EOL> gmean_sqrt = _map_eigenvalues ( np . sqrt , <EOL> conn_measure . mean_ ) <EOL> assert ( is_spd ( gmean_sqrt , decimal = <NUM_LIT:7> ) ) <EOL> assert ( is_spd ( conn_measure . whitening_ , decimal = <NUM_LIT:7> ) ) <EOL> assert_array_almost_equal ( conn_measure . whitening_ . dot ( <EOL> gmean_sqrt ) , np . eye ( n_features ) ) <EOL> assert_array_almost_equal ( gmean_sqrt . dot ( <EOL> _map_eigenvalues ( np . exp , cov_new ) ) . dot ( gmean_sqrt ) , <EOL> covs [ k ] ) <EOL> elif kind == "<STR_LIT>" : <EOL> assert ( is_spd ( cov_new , decimal = <NUM_LIT:7> ) ) <EOL> assert_array_almost_equal ( cov_new . dot ( covs [ k ] ) , <EOL> np . eye ( n_features ) ) <EOL> elif kind == "<STR_LIT>" : <EOL> assert ( is_spd ( cov_new , decimal = <NUM_LIT:7> ) ) <EOL> d = np . sqrt ( np . diag ( np . diag ( covs [ k ] ) ) ) <EOL> if cov_estimator == EmpiricalCovariance ( ) : <EOL> assert_array_almost_equal ( d . dot ( cov_new ) . dot ( d ) , <EOL> covs [ k ] ) <EOL> assert_array_almost_equal ( np . diag ( cov_new ) , <EOL> np . ones ( ( n_features ) ) ) <EOL> elif kind == "<STR_LIT>" : <EOL> prec = linalg . inv ( covs [ k ] ) <EOL> d = np . sqrt ( np . diag ( np . diag ( prec ) ) ) <EOL> assert_array_almost_equal ( d . dot ( cov_new ) . dot ( d ) , - prec + <EOL> <NUM_LIT:2> * np . diag ( np . diag ( prec ) ) ) </s>
