<s> from copy import deepcopy <EOL> import numpy as np <EOL> from scipy import linalg <EOL> from . . io . pick import pick_channels <EOL> from . . utils import logger , verbose <EOL> from . . forward import convert_forward_solution <EOL> from . . evoked import EvokedArray <EOL> from . . source_estimate import SourceEstimate <EOL> from . inverse import _subject_from_inverse <EOL> from . import apply_inverse <EOL> def _prepare_info ( inverse_operator ) : <EOL> """<STR_LIT>""" <EOL> info = deepcopy ( inverse_operator [ '<STR_LIT:info>' ] ) <EOL> info [ '<STR_LIT>' ] = <NUM_LIT> <EOL> info [ '<STR_LIT>' ] = inverse_operator [ '<STR_LIT>' ] <EOL> return info <EOL> def _pick_leadfield ( leadfield , forward , ch_names ) : <EOL> """<STR_LIT>""" <EOL> picks_fwd = pick_channels ( forward [ '<STR_LIT>' ] [ '<STR_LIT>' ] , ch_names ) <EOL> return leadfield [ picks_fwd ] <EOL> @ verbose <EOL> def point_spread_function ( inverse_operator , forward , labels , method = '<STR_LIT>' , <EOL> lambda2 = <NUM_LIT:1> / <NUM_LIT> , pick_ori = None , mode = '<STR_LIT>' , <EOL> n_svd_comp = <NUM_LIT:1> , verbose = None ) : <EOL> """<STR_LIT>""" <EOL> mode = mode . lower ( ) <EOL> if mode not in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] : <EOL> raise ValueError ( "<STR_LIT>" <EOL> % mode ) <EOL> logger . info ( "<STR_LIT>" % len ( labels ) ) <EOL> forward = convert_forward_solution ( forward , force_fixed = False , <EOL> surf_ori = True ) <EOL> info = _prepare_info ( inverse_operator ) <EOL> leadfield = _pick_leadfield ( forward [ '<STR_LIT>' ] [ '<STR_LIT:data>' ] [ : , <NUM_LIT:2> : : <NUM_LIT:3> ] , forward , <EOL> info [ '<STR_LIT>' ] ) <EOL> label_psf_summary = [ ] <EOL> label_singvals = [ ] <EOL> for ll in labels : <EOL> logger . info ( ll ) <EOL> if ll . hemi == '<STR_LIT>' : <EOL> offset = forward [ '<STR_LIT:src>' ] [ <NUM_LIT:0> ] [ '<STR_LIT>' ] . shape [ <NUM_LIT:0> ] <EOL> this_hemi = <NUM_LIT:1> <EOL> elif ll . hemi == '<STR_LIT>' : <EOL> offset = <NUM_LIT:0> <EOL> this_hemi = <NUM_LIT:0> <EOL> idx = np . intersect1d ( ll . vertices , forward [ '<STR_LIT:src>' ] [ this_hemi ] [ '<STR_LIT>' ] ) <EOL> fwd_idx = np . searchsorted ( forward [ '<STR_LIT:src>' ] [ this_hemi ] [ '<STR_LIT>' ] , idx ) <EOL> sub_leadfield = leadfield [ : , fwd_idx + offset ] <EOL> if mode == '<STR_LIT>' : <EOL> logger . info ( "<STR_LIT>" ) <EOL> this_label_psf_summary = sub_leadfield . sum ( axis = <NUM_LIT:1> ) [ np . newaxis , : ] <EOL> elif mode == '<STR_LIT>' : <EOL> logger . info ( "<STR_LIT>" ) <EOL> this_label_psf_summary = sub_leadfield . mean ( axis = <NUM_LIT:1> ) [ np . newaxis , : ] <EOL> elif mode == '<STR_LIT>' : <EOL> logger . info ( "<STR_LIT>" <EOL> % n_svd_comp ) <EOL> u_svd , s_svd , _ = linalg . svd ( sub_leadfield , <EOL> full_matrices = False , <EOL> compute_uv = True ) <EOL> label_singvals . append ( s_svd ) <EOL> logger . info ( "<STR_LIT>" % s_svd [ <NUM_LIT:0> : <NUM_LIT:5> ] ) <EOL> logger . info ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> my_comps = s_svd [ : n_svd_comp ] <EOL> comp_var = ( <NUM_LIT> * np . sum ( my_comps * my_comps ) / <EOL> np . sum ( s_svd * s_svd ) ) <EOL> logger . info ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( n_svd_comp , comp_var ) ) <EOL> this_label_psf_summary = ( u_svd [ : , : n_svd_comp ] * <EOL> s_svd [ : n_svd_comp ] [ np . newaxis , : ] ) <EOL> this_label_psf_summary = this_label_psf_summary . T <EOL> label_psf_summary . append ( this_label_psf_summary ) <EOL> label_psf_summary = np . concatenate ( label_psf_summary , axis = <NUM_LIT:0> ) <EOL> label_psf_summary = np . r_ [ label_psf_summary , <EOL> label_psf_summary . sum ( axis = <NUM_LIT:0> ) [ np . newaxis , : ] ] . T <EOL> evoked_fwd = EvokedArray ( label_psf_summary , info = info , tmin = <NUM_LIT:0.> ) <EOL> logger . info ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( method , lambda2 ) ) <EOL> stc_psf = apply_inverse ( evoked_fwd , inverse_operator , lambda2 , <EOL> method = method , pick_ori = pick_ori ) <EOL> return stc_psf , evoked_fwd <EOL> def _get_matrix_from_inverse_operator ( inverse_operator , forward , labels = None , <EOL> method = '<STR_LIT>' , lambda2 = <NUM_LIT:1.> / <NUM_LIT> , <EOL> mode = '<STR_LIT>' , n_svd_comp = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> mode = mode . lower ( ) <EOL> if not forward [ '<STR_LIT>' ] : <EOL> raise RuntimeError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> if not ( forward [ '<STR_LIT>' ] == <NUM_LIT:1> ) : <EOL> raise RuntimeError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> if labels : <EOL> logger . info ( "<STR_LIT>" % len ( labels ) ) <EOL> else : <EOL> logger . info ( "<STR_LIT>" ) <EOL> info = _prepare_info ( inverse_operator ) <EOL> id_mat = np . eye ( len ( info [ '<STR_LIT>' ] ) ) <EOL> ev_id = EvokedArray ( id_mat , info = info , tmin = <NUM_LIT:0.> ) <EOL> snr = <NUM_LIT> <EOL> lambda2 = <NUM_LIT:1.0> / snr ** <NUM_LIT:2> <EOL> invmat_mat_op = apply_inverse ( ev_id , inverse_operator , lambda2 = lambda2 , <EOL> method = method , pick_ori = '<STR_LIT>' ) <EOL> logger . info ( "<STR_LIT>" % str ( invmat_mat_op . shape ) ) <EOL> invmat_mat = invmat_mat_op . data <EOL> invmat_summary = [ ] <EOL> label_singvals = [ ] <EOL> if labels : <EOL> for ll in labels : <EOL> if ll . hemi == '<STR_LIT>' : <EOL> offset = forward [ '<STR_LIT:src>' ] [ <NUM_LIT:0> ] [ '<STR_LIT>' ] . shape [ <NUM_LIT:0> ] <EOL> this_hemi = <NUM_LIT:1> <EOL> elif ll . hemi == '<STR_LIT>' : <EOL> offset = <NUM_LIT:0> <EOL> this_hemi = <NUM_LIT:0> <EOL> else : <EOL> raise RuntimeError ( "<STR_LIT>" ) <EOL> idx = np . intersect1d ( ll . vertices , <EOL> forward [ '<STR_LIT:src>' ] [ this_hemi ] [ '<STR_LIT>' ] ) <EOL> fwd_idx = np . searchsorted ( forward [ '<STR_LIT:src>' ] [ this_hemi ] [ '<STR_LIT>' ] , idx ) <EOL> invmat_lbl = invmat_mat [ fwd_idx + offset , : ] <EOL> if mode == '<STR_LIT>' : <EOL> logger . info ( "<STR_LIT>" ) <EOL> this_invmat_summary = invmat_lbl . sum ( axis = <NUM_LIT:0> ) <EOL> this_invmat_summary = np . vstack ( this_invmat_summary ) . T <EOL> elif mode == '<STR_LIT>' : <EOL> logger . info ( "<STR_LIT>" ) <EOL> this_invmat_summary = invmat_lbl . mean ( axis = <NUM_LIT:0> ) <EOL> this_invmat_summary = np . vstack ( this_invmat_summary ) . T <EOL> elif mode == '<STR_LIT>' : <EOL> logger . info ( "<STR_LIT>" <EOL> "<STR_LIT>" % n_svd_comp ) <EOL> u_svd , s_svd , _ = linalg . svd ( invmat_lbl . T , <EOL> full_matrices = False , <EOL> compute_uv = True ) <EOL> label_singvals . append ( s_svd ) <EOL> logger . info ( "<STR_LIT>" % s_svd [ : <NUM_LIT:5> ] ) <EOL> logger . info ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> my_comps = s_svd [ : n_svd_comp ] <EOL> comp_var = ( ( <NUM_LIT:100> * np . sum ( my_comps * my_comps ) ) / <EOL> np . sum ( s_svd * s_svd ) ) <EOL> logger . info ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( n_svd_comp , comp_var ) ) <EOL> this_invmat_summary = ( u_svd [ : , : n_svd_comp ] . T * <EOL> s_svd [ : n_svd_comp ] [ : , np . newaxis ] ) <EOL> invmat_summary . append ( this_invmat_summary ) <EOL> invmat = np . concatenate ( invmat_summary , axis = <NUM_LIT:0> ) <EOL> else : <EOL> invmat = invmat_mat <EOL> return invmat , label_singvals <EOL> @ verbose <EOL> def cross_talk_function ( inverse_operator , forward , labels , <EOL> method = '<STR_LIT>' , lambda2 = <NUM_LIT:1> / <NUM_LIT> , signed = False , <EOL> mode = '<STR_LIT>' , n_svd_comp = <NUM_LIT:1> , verbose = None ) : <EOL> """<STR_LIT>""" <EOL> forward = convert_forward_solution ( forward , force_fixed = True , <EOL> surf_ori = True ) <EOL> out = _get_matrix_from_inverse_operator ( inverse_operator , forward , <EOL> labels = labels , method = method , <EOL> lambda2 = lambda2 , mode = mode , <EOL> n_svd_comp = n_svd_comp ) <EOL> invmat , label_singvals = out <EOL> leadfield = _pick_leadfield ( forward [ '<STR_LIT>' ] [ '<STR_LIT:data>' ] , forward , <EOL> inverse_operator [ '<STR_LIT:info>' ] [ '<STR_LIT>' ] ) <EOL> ctfs = np . dot ( invmat , leadfield ) <EOL> ctfs = np . vstack ( ( ctfs , ctfs . sum ( axis = <NUM_LIT:0> ) ) ) <EOL> if not signed : <EOL> ctfs = np . abs ( ctfs , out = ctfs ) <EOL> vertno = [ ss [ '<STR_LIT>' ] for ss in inverse_operator [ '<STR_LIT:src>' ] ] <EOL> stc_ctf = SourceEstimate ( ctfs . T , vertno , tmin = <NUM_LIT:0.> , tstep = <NUM_LIT:1.> ) <EOL> stc_ctf . subject = _subject_from_inverse ( inverse_operator ) <EOL> return stc_ctf </s>
