[Neurodebian-devel] Some meica.py fixes
Yaroslav Halchenko
debian at onerussian.com
Thu Nov 29 16:46:49 UTC 2012
Hi Prantik,
Please find attached a little patch up_meica_fixes which I
initiated while fixing an initial crash while trying to run:
meica.py -e 35 -d '06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_[1].nii.gz'
due to missing additional '%s' to match number of parameters provided
for string interpolation.
I also added few uses of os.path.join instead of "manual" %s/%s which
should be preferable in general, but feel free to skip those.
I started to look into those while trying to fix for a crash if I point
to the file laying in upstairs directory, i.e. in above using
'../06mar...' -- it would be nice if it was possible to point to files
at arbitrary location.
at the end of running above command it failed nevertheless with:
...
++ 3dTstat: AFNI version=AFNI_2011_12_21_1014 (Nov 28 2012) [64-bit]
++ Authored by: KR Hammett & RW Cox
++ Output dataset ./e1_std.nii.gz
-- ME-PCA/ME-ICA Component for ME-ICA v2.0 --
++ Loading Data
Traceback (most recent call last):
File "/home/yoh/deb/gits/pkg-exppsy/afni/pkundu/meica.libs/tedana.py", line 723, in <module>
catim = nib.load(options.data)
File "/usr/lib/pymodules/python2.7/nibabel/loadsave.py", line 39, in load
filename)
nibabel.spatialimages.ImageFileError: Cannot work out file type of "./06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3__in1"
and here is the content of directories at this stage:
$> ls
06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_1.nii.gz@ meica.06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_1/
_meica_06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_1.sh
$> ls meica.06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_1/
06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_1.nii.gz _eBvrmask.nii.gz eBvrmask.nii.gz
06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_1_vrA.1D _meica_06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_1.sh gms.1D
06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3_1_vrA.nii.gz e1_in.nii.gz motion.1D
06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3__vrmat.aff12.1D e1_mean.nii.gz stage
06mar11km_WIP_fMRI_SSh_3mm_sense2_sl35_SENSE_3__wmat.aff12.1D e1_std.nii.gz
_eBmask.nii.gz eBmask.nii.gz
so I just decided to check with you for further fixes needed to make
meica.py work.
I wondered -- may be it would be feasible to create a little "smoke
test" which would
1. generate some smallish 'fake' 4d sequence(s), e.g. using nibabel
so it would not be needed to drag some bulky 4d set for testing
2. run the script in few typical use-case scenarios
3. do basic tests of correct output, e.g. at least correct resolutions/sizes etc
?
exercising such a test would provide some guarantee that the script is
usable thus significantly reducing whining from users (like me ;) )
P.S. I also attached up_meica_tedana_upstream_mdp_api to which you would
need happen you decide to update MDP to the current state of its master
with your functionality adopted -- it is already available from
NeuroDebian.
Cheers
On Tue, 20 Nov 2012, Yaroslav Halchenko wrote:
> Thanks guys for the discussion. Although I have not touched ICAs for a
> while, having two separate limits for the coarse/fine optimizations
> makes sense in general. And it could be achieved without affecting
> existing users with only 1 small change -- not relying only on
> comparison of g and g_fine but rather on explicit setting for the
> complementary (coarser) limit. Since existing users had no clue about
> Initial and Secondary, and messages "Secondary convergence..." without
> ever seeing "Primary..." one, I suggest to make "primary_limit" a
> "coarse_limit" and set it to None, by default disabling this
> functionaliry. Then whoever needs to enable this coarse_limit (e.g. for
> meica) could do that.
> I have pushed original patch + my changes on top resulting in listed
> patch to:
> https://github.com/yarikoptic/mdp-toolkit
> and sent a PR
> I then have altered fastica tests to sweep also through having some
> coarse_limit set -- nothing crashes, although now it nearly doubles test run
> time, so you might like to improve on that ? ;) and I haven't tested with
> python3 yet.
> If this patch, in this or another shape gets accepted upstream
> in MDP then we could update MDP's debian package and use system-wide MDP for
> meica to avoid maintaining multiple copies/clones.
> Cheers,
> $> git diff origin/master..
> diff --git a/mdp/nodes/ica_nodes.py b/mdp/nodes/ica_nodes.py
> index 6ba1a0a..31f1657 100644
> --- a/mdp/nodes/ica_nodes.py
> +++ b/mdp/nodes/ica_nodes.py
> @@ -320,13 +320,14 @@ class FastICANode(ICANode):
> - 25.5.2005 now independent from scipy. Requires Numeric or numarray
> - 26.6.2006 converted to numpy
> - 14.9.2007 updated to Matlab version 2.5
> + - 26.6.2012 added ability to run two stages of optimization [PK]
> """
> def __init__(self, approach = 'defl', g = 'pow3', guess = None,
> fine_g = 'pow3', mu = 1,
> sample_size = 1, fine_tanh = 1, fine_gaus = 1,
> - max_it = 1000, max_it_fine = 100,
> - failures = 5, limit = 0.001, verbose = False,
> + max_it = 5000, max_it_fine = 100,
> + failures = 5, coarse_limit=None, limit = 0.001, verbose = False,
> whitened = False, white_comp = None, white_parm = None,
> input_dim = None, dtype=None):
> """
> @@ -346,6 +347,12 @@ def __init__(self, approach = 'defl', g = 'pow3', guess = None,
> It is passed directly to the WhiteningNode constructor.
> Ex: white_parm = { 'svd' : True }
> + coarse_limit -- initial convergence threshold, to switch to
> + fine_g function (i.e. linear to non-linear)
> + even before reaching the limit and final
> + tuning. Set it to a value higher than limit to
> + be in effect. PK 26-6-12.
> +
> limit -- convergence threshold.
> Specific for FastICA:
> @@ -417,6 +424,7 @@ def __init__(self, approach = 'defl', g = 'pow3', guess = None,
> self.fine_gaus = fine_gaus
> self.max_it = max_it
> self.max_it_fine = max_it_fine
> + self.coarse_limit = coarse_limit
> self.failures = failures
> self.guess = guess
> @@ -458,6 +466,7 @@ def core(self, data):
> guess = mult(guess, self.white.get_recmatrix(transposed=1))
> limit = self.limit
> + coarse_limit = self.coarse_limit
> max_it = self.max_it
> max_it_fine = self.max_it_fine
> failures = self.failures
> @@ -501,6 +510,7 @@ def core(self, data):
> used_g = gOrig
> stroke = 0
> fine_tuned = False
> + coarse_limit_reached = False
> lng = False
> # SYMMETRIC APPROACH
> @@ -529,6 +539,15 @@ def core(self, data):
> v2 = 1.-abs((mult(Q.T, QOldF)).diagonal()).min(axis=0)
> convergence_fine.append(v2)
> + if self.g != self.fine_g \
> + and coarse_limit is not None \
> + and convergence[round] < coarse_limit \
> + and not coarse_limit_reached:
> + if verbose:
> + print 'Coarse convergence, switching to fine cost...'
> + used_g = gFine
> + coarse_limit_reached = True
> +
> if convergence[round] < limit:
> if fine_tuning and (not fine_tuned):
> if verbose:
> @@ -569,7 +588,7 @@ def core(self, data):
> # Show the progress...
> if verbose:
> msg = ('Step no. %d,'
> - ' convergence: %.3f' % (round+1,convergence[round]))
> + ' convergence: %.7f' % (round+1,convergence[round]))
> print msg
> On Thu, 15 Nov 2012, Kundu, Prantik (NIH/NIMH) [F] wrote:
> > Hello All
> > As you have guessed, the two-stage ica_nodes.py in AFNI is a custom modification, intended to make possible one contrast to bootstrap another, for example enabling the use of tanh contrast in 'somewhat' noisy data by bootstrapping convergence with the more robust pow3 contrast. If primary and secondary are both given as the same contrast function, then the secondary step is skipped and the standard MDP algo continues. This modification is not in the MATLAB reference, and is something suited to the peculiarities of fMRI data. We know it improves sensitivity to spontaneous BOLD fluctuations in the course of denoising based on our NMR signal decay measures of BOLD signal quality, but otherwise it would be hard to tell the difference. Perhaps Tiziano would care to run this implementation on some reference data before deciding to integrate the patch into the mainline MDP. Importantly, the MDP as distributed in AFNI is isolated to the runtime directory of the code its distributed with (meica.libs), and is not made globally accessible.
> > -Prantik
> > ________________________________________
> > From: Tiziano Zito [tiziano.zito at bccn-berlin.de]
> > Sent: Thursday, November 15, 2012 3:52 PM
> > To: NeuroDebian Development; Kundu, Prantik (NIH/NIMH) [F]
> > Cc: mdp-dev
> > Subject: Re: MDP patches in AFNI
> > Hi,
> > I had a quick look at the AFNI localized copy of MDP. It seems they
> > imported a fairly recent git snapshot. If the ica_nodes.py
> > patch would get integrated in MDP, I think they could use the
> > released MDP 3.3 version without any incompatibility.
> > Is this 2 stage convergence switching also present in the original
> > MATLAB reference version of FastICA? I am far away from the ICA
> > world since a long time now, so I may have missed some important
> > developments lately. If this approach is indeed novel and introduced
> > by PK, wouldn't be better to allow the user to choose which
> > approach to use -- standard or 2-stage -- at instantiation time?
> > This would make MDP users happy, because their results are not going
> > to look different after upgrading to the next MDP release, and AFNI
> > could still use exactly the same algorithm their are using now...
> > What do you think?
> > Ciao,
> > Tiziano
> > PS: Ccing the mdp-dev mailing list.
> > On Thu 15 Nov, 15:08, Yaroslav Halchenko wrote:
> > > Hi Tiziano,
> > > NB Ccing PK -- original author of the patch
> > > As you might (or not) know, AFNI now carries a copy of MDP (of state
> > > somewhere before actual 3.3 release as far as I see) with slight
> > > patching of ica_nodes.py to have 2 stages of convergence with switching
> > > to gFine later in the loop (PK, please correct me if I am wrong)
> > > I am attaching patch for your consideration to get absorbed into MDP
> > > since we would like to avoid maintaining two copies of it ;)
> > > if needed (may be there is more patches, PK?) you can find
> > > complete AFNI source with this mdp copy at
> > > http://git.debian.org/?p=pkg-exppsy/afni.git
> > > branch upstream.
> > > --
> > > Yaroslav O. Halchenko
> > > Postdoctoral Fellow, Department of Psychological and Brain Sciences
> > > Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755
> > > Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419
> > > WWW: http://www.linkedin.com/in/yarik
> > > --- pkundu/meica.libs/mdp/nodes/ica_nodes_old.py 2012-11-15 14:33:53.000000000 -0500
> > > +++ pkundu/meica.libs/mdp/nodes/ica_nodes.py 2012-11-15 14:33:53.000000000 -0500
> > > @@ -1,3 +1,4 @@
> > > +
> > > __docformat__ = "restructuredtext en"
> > > import math
> > > @@ -325,8 +326,8 @@
> > > def __init__(self, approach = 'defl', g = 'pow3', guess = None,
> > > fine_g = 'pow3', mu = 1,
> > > sample_size = 1, fine_tanh = 1, fine_gaus = 1,
> > > - max_it = 1000, max_it_fine = 100,
> > > - failures = 5, limit = 0.001, verbose = False,
> > > + max_it = 5000, max_it_fine = 100,
> > > + failures = 5, primary_limit=0.01, limit = 0.001, verbose = False,
> > > whitened = False, white_comp = None, white_parm = None,
> > > input_dim = None, dtype=None):
> > > """
> > > @@ -346,7 +347,9 @@
> > > It is passed directly to the WhiteningNode constructor.
> > > Ex: white_parm = { 'svd' : True }
> > > - limit -- convergence threshold.
> > > + limit -- final convergence threshold.
> > > +
> > > + primary_limit -- initial convergence threshold, to switch to fine function, (i.e. linear to non-linear). PK 26-6-12.
> > > Specific for FastICA:
> > > @@ -417,6 +420,7 @@
> > > self.fine_gaus = fine_gaus
> > > self.max_it = max_it
> > > self.max_it_fine = max_it_fine
> > > + self.primary_limit = primary_limit
> > > self.failures = failures
> > > self.guess = guess
> > > @@ -458,6 +462,7 @@
> > > guess = mult(guess, self.white.get_recmatrix(transposed=1))
> > > limit = self.limit
> > > + primary_limit = self.primary_limit
> > > max_it = self.max_it
> > > max_it_fine = self.max_it_fine
> > > failures = self.failures
> > > @@ -501,6 +506,7 @@
> > > used_g = gOrig
> > > stroke = 0
> > > fine_tuned = False
> > > + in_secondary = False
> > > lng = False
> > > # SYMMETRIC APPROACH
> > > @@ -529,10 +535,16 @@
> > > v2 = 1.-abs((mult(Q.T, QOldF)).diagonal()).min(axis=0)
> > > convergence_fine.append(v2)
> > > + if self.g!=self.fine_g and convergence[round] < primary_limit and not in_secondary:
> > > + if verbose:
> > > + print 'Primary convergence, switching to fine cost...'
> > > + used_g = gFine
> > > + in_secondary = True
> > > +
> > > if convergence[round] < limit:
> > > if fine_tuning and (not fine_tuned):
> > > if verbose:
> > > - print 'Initial convergence, fine-tuning...'
> > > + print 'Secondary convergence, fine-tuning...'
> > > fine_tuned = True
> > > used_g = gFine
> > > mu = muK * self.mu
> > > @@ -569,7 +581,7 @@
> > > # Show the progress...
> > > if verbose:
> > > msg = ('Step no. %d,'
> > > - ' convergence: %.3f' % (round+1,convergence[round]))
> > > + ' convergence: %.7f' % (round+1,convergence[round]))
> > > print msg
> > _______________________________________________
> > Neurodebian-devel mailing list
> > Neurodebian-devel at lists.alioth.debian.org
> > http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/neurodebian-devel
--
Yaroslav O. Halchenko
Postdoctoral Fellow, Department of Psychological and Brain Sciences
Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755
Phone: +1 (603) 646-9834 Fax: +1 (603) 646-1419
WWW: http://www.linkedin.com/in/yarik
-------------- next part --------------
From: Yaroslav Halchenko <debian at onerussian.com>
Subject: Initiated by fixing crash while providing a single 4D volume
since it requires 1 more formatcard (%s) in the string for the interpolation
also complemented with use of os.path.merge as exemplars to write OS-agnostic
path handling
Origin: NeuroDebian
Last-Update: 2012-11-29
--- a/pkundu/meica.py
+++ b/pkundu/meica.py
@@ -283,7 +283,7 @@ if options.anat!='':
sl.append("3dcalc -a eBmask%s -b eBmask_bias%s -expr 'a/b' -prefix eBbase%s" % ( osf, osf, osf))
weightline = ' -lpc -weight eBmask_pve_0.nii.gz -base eBbase%s ' % (osf)
else: weightline = ' -lpc+ZZ -automask -autoweight -base _eBmask%s ' % (osf)
- sl.append("cp %s/%s* ." % (startdir,nsmprage))
+ sl.append("cp %s* ." % os.path.merge(startdir,nsmprage))
abmprage = nsmprage
if options.tlrc:
sl.append("afnibinloc=`which 3dSkullStrip`")
@@ -294,7 +294,7 @@ if options.anat!='':
sl.append("3dcopy %s %s" % (atnsmprage,dsprefix(atnsmprage)))
sl.append("3drefit -view orig %s+tlrc " % dsprefix(atnsmprage) )
sl.append("cp %s %s" % (atnsmprage,startdir))
- sl.append("gzip -f %s/%s" % (startdir,atnsmprage))
+ sl.append("gzip -f %s" % os.path.merge(startdir,atnsmprage))
abmprage = atnsmprage
align_args=""
@@ -302,7 +302,7 @@ if options.anat!='':
elif oblique_mode: align_args = " -maxrot 10 -maxshf 10 "
else: align_args=" -maxrot 20 -maxshf 20 -parfix 7 1 -parang 9 0.83 1.0 "
if oblique_mode: alnsmprage = "./%s_ob.nii.gz" % (anatprefix)
- else: alnsmprage = "%s/%s" % (startdir,nsmprage)
+ else: alnsmprage = os.path.merge(startdir,nsmprage)
sl.append("3dAllineate -weight_frac 1.0 -VERB -warp aff %s -source_automask+10 -cmass -master SOURCE -source %s -prefix ./%s_al -1Dmatrix_save %s_al_mat %s" % (weightline,alnsmprage, anatprefix,anatprefix,align_args))
if options.tlrc: tlrc_opt = "%s::WARP_DATA -I" % (atnsmprage)
else: tlrc_opt = ""
@@ -366,7 +366,7 @@ sl.append("echo 2 > stage")
if len(ica_datasets)==1:
dsin = ''.join(ica_datasets)+trailing
ica_prefix = dsin
- ica_input="./%s_in%s" % (prefix,''.join(ica_datasets),trailing)
+ ica_input="./%s_in%s%s" % (prefix,''.join(ica_datasets),trailing)
ica_mask="eBvrmask.nii.gz"
else:
ica_input = "zcat_ffd.nii.gz"
-------------- next part --------------
From: Yaroslav Halchenko <debian at onerussian.com>
Subject: Use coarse_limit in FastICANode parameters, as patch was adopted upstream
Corresponding patch from meica's copy of MDP was adopted upstream in MDP
as of MDP-3.3-6-g7bbd889 :
https://github.com/mdp-toolkit/mdp-toolkit/commit/7bbd8896edf2de30d4c6376d1a37b02cc2a94023
Last-Update: 2012-11-28
--- a/pkundu/meica.libs/tedana.py
+++ b/pkundu/meica.libs/tedana.py
@@ -580,7 +580,7 @@ def tedica(dd,cost):
#Do ICA
climit = float("%s" % options.conv)
#icanode = mdp.nodes.FastICANode(white_comp=nc, white_parm={'svd':True},approach='symm', g=cost, fine_g=options.finalcost, limit=climit, verbose=True)
- icanode = mdp.nodes.FastICANode(white_comp=nc,approach='symm', g=cost, fine_g=options.finalcost, primary_limit=climit*100, limit=climit, verbose=True)
+ icanode = mdp.nodes.FastICANode(white_comp=nc,approach='symm', g=cost, fine_g=options.finalcost, coarse_limit=climit*100, limit=climit, verbose=True)
icanode.train(dd)
smaps = icanode.execute(dd)
mmix = icanode.get_recmatrix().T
More information about the Neurodebian-devel
mailing list