Thanks for that Gael - I do know nilearn but in this case I did a depth-first search on doing SVM on brain images and ended up here :) The size of 'coef' is (1,205739), the size of mask[mask].T is (205739,)
Same number of elements, different storage layout(?). Turned out that the numpy.ravel() function -pointed out by my colleague- can solve that! >>> wmap[mask]=np.ravel(coef) bw Alle Meije On 5 May 2017 at 07:59, Gael Varoquaux <gael.varoqu...@normalesup.org> wrote: > Hi Alle, > > I think that what has changed between 2014 and today is that the > coefficients coef are now a 2D array (number of hyperplanes x number of > features). In your case, the first direction is of length one, so you > could just do: > > coef = clf.coef_[0] > > and your script should work. > > The code of the Abraham paper cannot be updated, because papers are not > living objects. However, we are maintaining a package that encodes all > these patterns in higher-level construct: http://nilearn.github.io/ > It might be a good idea to use this package, as it is maintained and has > quality assurance. > > Best, > > Gaël > > On Thu, May 04, 2017 at 05:52:27PM +0200, Alle Meije Wink wrote: > > I have a script to classify MRI perfusion maps from healthy subjects and > > patients. For the file IO and the classifier I have started with the > example > > code in Abraham et al 2014 [https://arxiv.org/pdf/1412.3919.pdf]. > > > I use the same classifier as in the paper to produce a back-projected > map of > > classification weights, which I then want to 'unmask' like in the paper: > > > coef=clf.coef_ > > coef=featureselection.inverse_transform(coef) > > > and > > > map_name='weights_check.nii.gz' > > wmap=np.zeros(mask.shape, dtype=X.dtype) > > wmap[mask]=coef > > img=nb.Nifti1Image(wmap,np.eye(4)) > > img.to_filename(map_name) > > > But the line "wmap[mask]=coef" throws an error "ValueError: boolean > index array > > should have 1 dimension". I tried the example code from the paper and > that > > works. > > > Is the 'coef' array of back-projected SVM weights in some way different > than > > the masked input image? Or am I doing something else wrong? The error > suggests > > that the mask array is the problem. > > > The complete script is attached. > > > Many thanks for your help! > > Alle Meije > > > # -*- coding: utf-8 -*- > > """ > > classify MR images from controls (CON) and patients (MCI) > > > """ > > > import os > > import numpy as np > > import nibabel as nb > > import sklearn as sl > > > from sklearn.feature_selection import f_classif > > from sklearn import svm > > > def subj_lists( rootdir='/data/a.wink/MR', starts=['CON','MCI'] ): > > > # make and empty list for each patient group > > slists=[] > > for i in range(0,len(starts)): > > slists.append([]); > > > # add subjects to the group based on > > for root, dirs, files in os.walk(rootdir): > > for fname in files: > > for i in range(0,len(starts)): > > if fname.startswith(starts[i]): > > slists[i].append(os.path.join(root,fname)) > > > print('finished building subject lists') > > return slists > > > def mk_mask ( subj_lists=[] ): > > > # check whether mask can be loaded > > mname='mask_check.nii.gz' > > > if os.path.isfile(mname): > > > msk=nb.load(mname).get_data() > > > # or must be built (values >20% in the subject images and >20% grey > matter in the template) > > else: > > > num_im=0; > > for i in range(0,len(subj_lists)): > > for fname in subj_lists[i]: > > imdata=nb.load(fname).get_data() > > if 'msk' not in locals(): > > msk=np.absolute(imdata) > > else: > > msk=msk+np.absolute(imdata) > > num_im+=1 > > > msk/=num_im > > msk/=np.amax(msk) > > msk[msk<.2]=0 > > msk[msk>0]=1 > > > grey=nb.load('/usr/local/spm8/apriori/grey.nii').get_data() > > grey[grey<.2]=0 > > grey[grey>0]=1 > > > msk=msk*grey > > > img=nb.Nifti1Image(msk,np.eye(4)) > > img.to_filename(mname) > > > msk=msk.astype(bool) > > > print('finished building mask') > > return msk, mname > > > def load_images( subj_lists=[], mask=[] ): > > > # load all the images, build a matrix X or subjects (rows) * > inmask-voxels (columns) > > num_im=0; > > for i in range(0,len(subj_lists)): > > for fname in subj_lists[i]: > > imdata=nb.load(fname).get_data() > > num_im+=1 > > print '\r' + str(num_im) +'\r', > > if 'the_matrix' not in locals(): > > the_matrix=imdata[mask].T > > else: > > the_matrix=np.vstack((the_matrix,imdata[mask].T)) > > > print('finished building matrix X') > > return the_matrix > > > def main(): > > > # get the input data > > subjlists = subj_lists() > > y=np.concatenate( ([-1 for _ in range(len(subjlists[0]))], > > [ 1 for _ in range(len(subjlists[1]))]) ) > > mask, mas_name = mk_mask(subjlists) > > X = load_images(subjlists,mask) > > print "number of subjects, voxels: %d, %d" % X.shape > > > # select features > > featureselection=sl.feature_selection.SelectKBest(f_classif, k=8000) > > X_reduced=featureselection.fit_transform(X,y) > > > # classify > > clf=svm.SVC(kernel='linear') > > clf.fit(X_reduced,y) > > > # make discrimination map > > coef=clf.coef_ > > coef=featureselection.inverse_transform(coef) > > > map_name='weights_check.nii.gz' > > wmap=np.zeros(mask.shape, dtype=X.dtype) > > wmap[mask]=coef > > img=nb.Nifti1Image(wmap,np.eye(4)) > > img.to_filename(map_name) > > > if __name__ == "__main__": > > main() > > > > _______________________________________________ > > scikit-learn mailing list > > scikit-learn@python.org > > https://mail.python.org/mailman/listinfo/scikit-learn > > > -- > Gael Varoquaux > Researcher, INRIA Parietal > NeuroSpin/CEA Saclay , Bat 145, 91191 Gif-sur-Yvette France > Phone: ++ 33-1-69-08-79-68 > http://gael-varoquaux.info http://twitter.com/GaelVaroquaux > _______________________________________________ > scikit-learn mailing list > scikit-learn@python.org > https://mail.python.org/mailman/listinfo/scikit-learn >
_______________________________________________ scikit-learn mailing list scikit-learn@python.org https://mail.python.org/mailman/listinfo/scikit-learn
