<div dir="ltr"><span style="font-size:12.8px">Many thanks for all the feedback. I get the highest accuracy at dimension [15, 45, 65], not sure if everything is correct - logically in the code now. I am yet classifying whole brain, just as a check.</span><div><span style="font-size:12.8px">Also, the plotting seems completely off...and the output either non-existent or really weird. Is there any chance anyone could quickly run the code? Or would you have any other suggestions?</span><br></div><div><div style="font-size:12.8px">Best regards,</div><div style="font-size:12.8px">Pegah</div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px"><b>- Newest code below</b></div><div><b style="font-size:12.8px">- Data is here: </b><span style="font-size:12.8px"><b><a href="https://www.dropbox.com/sh/6qnrt2l2othc83g/AADBpc-eJSK5893Vrz_A8BS1a?dl=0">https://www.dropbox.com/sh/6qnrt2l2othc83g/AADBpc-eJSK5893Vrz_A8BS1a?dl=0</a></b></span></div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px">Code:</div><div style="font-size:12.8px"><br></div><div style="font-size:12.8px"><div>#</div><span class="gmail-im"><div><br></div><div>from glob import glob</div><div>import os</div><div>import numpy as np</div><div><br></div><div>from mvpa2.suite import *</div><div><br></div><div>%matplotlib inline</div><div><br></div><div><br></div><div># enable debug output for searchlight call</div><div>if __debug__:</div><div> debug.active += ["SLC"]</div><div><br></div><div><br></div></span><div># change working HERE!</div><div>os.chdir('mypath/S1')</div><span class="gmail-im"><div><br></div><div># use glob to get the filenames of .nii data into a list</div><div>nii_fns = glob('beta*.nii')</div><div><br></div><div><br></div></span><span class="gmail-im"><div># read data</div><div><br></div><div>labels = [</div><div> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, </div><div> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, </div><div> 7, 7, 7, 7, 7, 7,</div><div> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, </div><div> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, </div></span><div> 7, 7, 7, 7, 7, 7, 7, 7</div><div>]</div><span class="gmail-im"><div>grps = np.repeat([0, 1], 37, axis=0) # used for `chuncks`</div><div><br></div></span><div>sprefix_ = 'vxl'</div><div>sprefix_indices_key = '_'.join([sprefix_, 'indices'])</div><div>tprefix_ = 'tpref'</div><div>db = mvpa2.datasets.mri.fmri_datase<wbr>t(</div><div> nii_fns, targets=labels, chunks=grps, mask=None, sprefix=sprefix_, tprefix=tprefix_, add_fa=None</div><span class="gmail-im"><div>)</div><div><br></div><div><br></div><div># use only the samples of which labels are 1 or 2</div></span><span class="gmail-im"><div>db12 = db.select(sadict={'targets': [1,2]})</div><div><br></div><div># in-place z-score normalization</div><div>zscore(db12)</div><div><br></div><div><br></div><div># choose classifier</div><div>clf = LinearNuSVMC()</div><div><br></div><div># setup measure to be computed by Searchlight</div><div># cross-validated mean transfer using an N-fold dataset splitter</div><div>cv = CrossValidation(clf, NFoldPartitioner())</div><div><br></div><div># define searchlight methods</div></span><div>radius_ = 3</div><div>sl = sphere_searchlight(</div><div> cv, radius=radius_, space=sprefix_indices_key,</div><span class="gmail-im"><div> postproc=mean_sample()</div><div>)</div><div><br></div><div># stripping all attributes from the dataset that are not required for the searchlight analysis</div><div>db12s = db12.copy(</div><div> deep=False,</div><div> sa=['targets', 'chunks'],</div></span><div> fa=[sprefix_indices_key],</div><span class="gmail-im"><div> a=['mapper']</div><div>)</div><div><br></div><div># run searchlight</div><div>sl_map = sl(db12s)</div><div><br></div><div><br></div><div># toggle between error rate (searchlight output) and accuracy (for plotting)</div><div>sl_map.samples *= -1</div><div>sl_map.samples += 1</div><div><br></div><div># The result dataset is fully aware of the original dataspace. </div><div># Using this information we can map the 1D accuracy maps back into “brain-space” (using NIfTI image header information from the original input timeseries.</div><div>niftiresults = map2nifti(sl_map, imghdr=db12.a.imghdr)</div><div><br></div><div><br></div><div># find the 3-d coordinates of extrema of error rate or accuracy</div><div>extremum_i = np.argmax(sl_map.samples[0]) # max</div><div>extremum_i = np.argmin(sl_map.samples[0]) # min</div></span><div>coord = db12s.fa[sprefix_indices_key][<wbr>extremum_i]</div><span class="gmail-im"><div><br></div><div><br></div><div># plotting</div><div><br></div><div>plot_args = {</div><div> 'background' : 'Struct.nii',</div><div> 'background_mask' : 'brain.nii',</div></span><div> 'overlay_mask' : 'S1_mask.nii',</div><span class="gmail-im"><div> 'do_stretch_colors' : False,</div><div> 'cmap_bg' : 'gray',</div><div> 'cmap_overlay' : 'autumn', # YlOrRd_r # pl.cm.autumn</div><div> 'interactive' : cfg.getboolean('examples', 'interactive', True)</div><div>}</div><div><br></div><div>fig = pl.figure(figsize=(12, 4), facecolor='white')</div><div><br></div><div>subfig = plot_lightbox(</div><div> overlay=niftiresults,</div><div> vlim=(0.5, None), slices=range(23,31),</div><div> fig=fig, </div><div> background=Struct.ni<wbr>i',</div></span><div> background_mask='brain.nii',</div><span class="gmail-im"><div> overlay_mask='S1_mask.nii',</div><div> **plot_args</div><div>)</div><div><br></div><div>pl.title('Accuracy distribution for radius %i' % radius_)</div><div><br></div><div><br></div></span><div>#</div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sat, Sep 9, 2017 at 5:29 PM, Pegah Kassraian Fard <span dir="ltr"><<a href="mailto:pegahkf@gmail.com" target="_blank">pegahkf@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Many thanks for all the feedback. I get the highest accuracy at dimension [15, 45, 65], not sure if everything is correct - logically in the code now - the plotting seems completely off. Is there any chance anyone could quickly run the code? Or would you have any other suggestions?<div>Best regards,</div><div>Pegah</div><div><br></div><div><b>Attached:</b></div><div><b>- Newest code</b></div><div><b>- Output</b></div><div><b>- Plot</b></div><div><br></div><div><br></div><div>Code also here:</div><div><br></div><div><div>#</div><span class=""><div><br></div><div>from glob import glob</div><div>import os</div><div>import numpy as np</div><div><br></div><div>from mvpa2.suite import *</div><div><br></div><div>%matplotlib inline</div><div><br></div><div><br></div><div># enable debug output for searchlight call</div><div>if __debug__:</div><div> debug.active += ["SLC"]</div><div><br></div><div><br></div></span><div># change working directory to 'S1'</div><div>os.chdir('mypath/S1')</div><span class=""><div><br></div><div># use glob to get the filenames of .nii data into a list</div><div>nii_fns = glob('beta*.nii')</div><div><br></div><div><br></div></span><span class=""><div># read data</div><div><br></div><div>labels = [</div><div> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, </div><div> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, </div><div> 7, 7, 7, 7, 7, 7,</div><div> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, </div><div> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, </div></span><div> 7, 7, 7, 7, 7, 7, 7, 7</div><div>]</div><span class=""><div>grps = np.repeat([0, 1], 37, axis=0) # used for `chuncks`</div><div><br></div></span><div>sprefix_ = 'vxl'</div><div>sprefix_indices_key = '_'.join([sprefix_, 'indices'])</div><div>tprefix_ = 'tpref'</div><div>db = mvpa2.datasets.mri.fmri_datase<wbr>t(</div><div> nii_fns, targets=labels, chunks=grps, mask=None, sprefix=sprefix_, tprefix=tprefix_, add_fa=None</div><span class=""><div>)</div><div><br></div><div><br></div><div># use only the samples of which labels are 1 or 2</div></span><span class=""><div>db12 = db.select(sadict={'targets': [1,2]})</div><div><br></div><div># in-place z-score normalization</div><div>zscore(db12)</div><div><br></div><div><br></div><div># choose classifier</div><div>clf = LinearNuSVMC()</div><div><br></div><div># setup measure to be computed by Searchlight</div><div># cross-validated mean transfer using an N-fold dataset splitter</div><div>cv = CrossValidation(clf, NFoldPartitioner())</div><div><br></div><div># define searchlight methods</div></span><div>radius_ = 3</div><div>sl = sphere_searchlight(</div><div> cv, radius=radius_, space=sprefix_indices_key,</div><span class=""><div> postproc=mean_sample()</div><div>)</div><div><br></div><div># stripping all attributes from the dataset that are not required for the searchlight analysis</div><div>db12s = db12.copy(</div><div> deep=False,</div><div> sa=['targets', 'chunks'],</div></span><div> fa=[sprefix_indices_key],</div><span class=""><div> a=['mapper']</div><div>)</div><div><br></div><div># run searchlight</div><div>sl_map = sl(db12s)</div><div><br></div><div><br></div><div># toggle between error rate (searchlight output) and accuracy (for plotting)</div><div>sl_map.samples *= -1</div><div>sl_map.samples += 1</div><div><br></div><div># The result dataset is fully aware of the original dataspace. </div><div># Using this information we can map the 1D accuracy maps back into “brain-space” (using NIfTI image header information from the original input timeseries.</div><div>niftiresults = map2nifti(sl_map, imghdr=db12.a.imghdr)</div><div><br></div><div><br></div><div># find the 3-d coordinates of extrema of error rate or accuracy</div><div>extremum_i = np.argmax(sl_map.samples[0]) # max</div><div>extremum_i = np.argmin(sl_map.samples[0]) # min</div></span><div>coord = db12s.fa[sprefix_indices_key][<wbr>extremum_i]</div><span class=""><div><br></div><div><br></div><div># plotting</div><div><br></div><div>plot_args = {</div><div># 'background' : 'Subject2_Struct.nii',</div><div># 'background_mask' : 'brain.nii.gz',</div></span><div># 'overlay_mask' : 'S1.nii',</div><span class=""><div> 'do_stretch_colors' : False,</div><div> 'cmap_bg' : 'gray',</div><div> 'cmap_overlay' : 'autumn', # YlOrRd_r # pl.cm.autumn</div><div> 'interactive' : cfg.getboolean('examples', 'interactive', True)</div><div>}</div><div><br></div><div>fig = pl.figure(figsize=(12, 4), facecolor='white')</div><div><br></div><div>subfig = plot_lightbox(</div><div> overlay=niftiresults,</div><div> vlim=(0.5, None), slices=range(23,31),</div><div> fig=fig, </div><div> background='Subject2_Struct.ni<wbr>i',</div></span><div># background_mask='brain.nii',</div><span class=""><div> overlay_mask='S1_mask.nii',</div><div> **plot_args</div><div>)</div><div><br></div><div>pl.title('Accuracy distribution for radius %i' % radius_)</div><div><br></div><div><br></div><div><br></div><div><br></div></span><div>#</div></div><div><br></div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Sep 8, 2017 at 6:26 PM, Yaroslav Halchenko <span dir="ltr"><<a href="mailto:debian@onerussian.com" target="_blank">debian@onerussian.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span><br>
On Fri, 08 Sep 2017, Pegah Kassraian Fard wrote:<br>
<br>
<br>
<br>
<br>
</span><span>> from glob import glob<br>
> import os<br>
> import numpy as np<br>
<br>
> from mvpa2.suite import *<br>
<br>
> %matplotlib inline<br>
<br>
<br>
> # enable debug output for searchlight call<br>
> if __debug__:<br>
> debug.active += ["SLC"]<br>
<br>
<br>
> # change working directory to 'WB'<br>
> os.chdir('mypath/WB')<br>
<br>
> # use glob to get the filenames of .nii data into a list<br>
> nii_fns = glob('beta*.nii')<br>
<br>
</span>glob order iirc might not be guaranteed, so I would sort it to make sure<br>
<br>
nii_fns = sorted(glob('beta*.nii'))<br>
<span><br>
> # read data<br>
<br>
> labels = [<br>
> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,<br>
> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,<br>
> 7, 7, 7, 7, 7, 7, 7,<br>
> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,<br>
> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,<br>
> 7, 7, 7, 7, 7, 7, 7<br>
> ]<br>
<br>
</span>and make sure that above order of the volumes correspond to those labels<br>
<span><br>
> grps = np.repeat([0, 1], 37, axis=0) # used for `chuncks`<br>
<br>
> db = mvpa2.datasets.mri.fmri_datase<wbr>t(<br>
> nii_fns, targets=labels, chunks=grps, mask=None, sprefix='vxl', tprefix='tpref', add_fa=None<br>
> )<br>
<br>
> # use only the samples of which labels are 1 or 2<br>
> db12 = db[np.array([label in [1, 2] for label in labels], dtype='bool')]<br>
<br>
</span>more concise way<br>
<br>
db12 = db.select(sadict={'targets': [1,2]})<br>
<div><div class="m_-6878619627213120064h5"><br>
> # in-place z-score normalization<br>
> zscore(db12)<br>
<br>
> # choose classifier<br>
> clf = LinearNuSVMC()<br>
<br>
> # setup measure to be computed by Searchlight<br>
> # cross-validated mean transfer using an N-fold dataset splitter<br>
> cv = CrossValidation(clf, NFoldPartitioner())<br>
<br>
> # define searchlight methods<br>
> radius_ = 1<br>
> sl = sphere_searchlight(<br>
> cv, radius=radius_, space='vxl_indices',<br>
> postproc=mean_sample()<br>
> )<br>
<br>
> # stripping all attributes from the dataset that are not required for the searchlight analysis<br>
> db12s = db12.copy(<br>
> deep=False,<br>
> sa=['targets', 'chunks'],<br>
> fa=['vxl_indices'],<br>
> a=['mapper']<br>
> )<br>
<br>
> # run searchlight<br>
> sl_map = sl(db12s)<br>
<br>
> # toggle between error rate (searchlight output) and accuracy (for plotting)<br>
> sl_map.samples *= -1<br>
> sl_map.samples += 1<br>
<br>
> # The result dataset is fully aware of the original dataspace.<br>
> # Using this information we can map the 1D accuracy maps back into “brain-space” (using NIfTI image header information from the original input timeseries.<br>
> niftiresults = map2nifti(sl_map, imghdr=db12.a.imghdr)<br>
<br>
> # find the 3-d coordinates of extrema of error rate or accuracy<br>
> extremum_i = np.argmax(sl_map.samples[0]) # max<br>
> extremum_i = np.argmin(sl_map.samples[0]) # min<br>
<br>
</div></div>note that above you override the same extremum_i with 'min' (minimal accuracy ;)) !!!<br>
<span><br>
> coord = db12s.fa[list(db.fa.keys())[0]<wbr>][extremum_i]<br>
<br>
</span>.keys() order is also not guaranteed... why not to specify<br>
'vxl_indices'?<br>
<span class="m_-6878619627213120064im m_-6878619627213120064HOEnZb"><br>
> # plotting<br>
> plot_args = {<br>
> # 'background' : 'Subject2_Struct.nii',<br>
> # 'background_mask' : 'brain.nii.gz',<br>
> # 'overlay_mask' : 'S1_mask.nii',<br>
> 'do_stretch_colors' : False,<br>
> 'cmap_bg' : 'gray',<br>
> 'cmap_overlay' : 'autumn', # YlOrRd_r # pl.cm.autumn<br>
> 'interactive' : cfg.getboolean('examples', 'interactive', True)<br>
> }<br>
> fig = pl.figure(figsize=(12, 4), facecolor='white')<br>
> subfig = plot_lightbox(<br>
> overlay=niftiresults,<br>
> vlim=(0.5, None), slices=range(23,31),<br>
> fig=fig,<br>
> background='Subject2_Struct.n<wbr>ii',<br>
> # background_mask='brain.nii.<wbr>gz',<br>
> # overlay_mask='S1_mask.nii',<br>
> **plot_args<br>
> )<br>
> pl.title('Accuracy distribution for radius %i' % radius_)<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
</span><div class="m_-6878619627213120064HOEnZb"><div class="m_-6878619627213120064h5">--<br>
Yaroslav O. Halchenko<br>
Center for Open Neuroscience <a href="http://centerforopenneuroscience.org" rel="noreferrer" target="_blank">http://centerforopenneuroscie<wbr>nce.org</a><br>
Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755<br>
Phone: <a href="tel:%2B1%20%28603%29%20646-9834" value="+16036469834" target="_blank">+1 (603) 646-9834</a> Fax: <a href="tel:%2B1%20%28603%29%20646-1419" value="+16036461419" target="_blank">+1 (603) 646-1419</a><br>
WWW: <a href="http://www.linkedin.com/in/yarik" rel="noreferrer" target="_blank">http://www.linkedin.com/in/ya<wbr>rik</a><br>
<br>
______________________________<wbr>_________________<br>
Pkg-ExpPsy-PyMVPA mailing list<br>
<a href="mailto:Pkg-ExpPsy-PyMVPA@lists.alioth.debian.org" target="_blank">Pkg-ExpPsy-PyMVPA@lists.alioth<wbr>.debian.org</a><br>
<a href="http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-exppsy-pymvpa" rel="noreferrer" target="_blank">http://lists.alioth.debian.org<wbr>/cgi-bin/mailman/listinfo/pkg-<wbr>exppsy-pymvpa</a></div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>