<commit_before>import gc, os, fpformat
from neuroimaging import traits

import numpy as N
import numpy.linalg as L
import numpy.random as R
import scipy.ndimage
from scipy.sandbox.models.utils import monotone_fn_inverter, rank 
from scipy.sandbox.models.regression import OLSModel, ARModel

from neuroimaging.modalities.fmri import fMRIImage
from neuroimaging.modalities.fmri.fmristat.delay import DelayContrast, DelayContrastOutput
from neuroimaging.modalities.fmri.regression import AROutput, TContrastOutput, \
  FContrastOutput, ResidOutput
from neuroimaging.algorithms import kernel_smooth
from neuroimaging.algorithms.fwhm import fastFWHM
from neuroimaging.algorithms.utils import fwhm2sigma
from neuroimaging.algorithms.statistics.regression import LinearModelIterator


# Imports related to pylab

from neuroimaging.defines import pylab_def
PYLAB_DEF, pylab = pylab_def()
if PYLAB_DEF:
    from neuroimaging.ui.visualization.multiplot import MultiPlot

class WholeBrainNormalize(traits.HasTraits):

    mask = traits.Any()

    def __init__(self, fmri_image, **keywords):
        traits.HasTraits.__init__(self, **keywords)
        if self.mask is not None:
            self._mask = self.mask.readall()
            self._mask.shape = N.product(self._mask.shape)
            
        self.n = fmri_image.grid.shape[0]
        self.avg = N.zeros((self.n,), N.float64)

        for i in range(self.n):
            d = fmri_image[slice(i,i+1)]
            if hasattr(self, '_mask'):
                d.shape = N.product(d.shape)
                d = N.compress(self._mask, d)
            self.avg[i] = d.mean()

    def __call__(self, fmri_data):
        out = N.zeros(fmri_data.shape, N.float64)
        for i in range(self.n):
            out[i] = fmri_data[i] * 100. / self.avg[i]
        return out

class fMRIStatOLS(LinearModelIterator):

    """
    OLS pass of fMRIstat.
    """
    
    normalize = traits.Any()
    
    formula = traits.Any()
    slicetimes = traits.Any()
    tshift = traits.Float(0.)

    fwhm_rho = traits.Float(6.)
    fwhm_data = traits.Float(6.)
    target_df = traits.Float(100.)
    nmax = traits.Int(200) # maximum number of rho values
    mask = traits.Any()
    path = traits.String('fmristat_run')
    resid = traits.false
    clobber = traits.false
    output_fwhm = traits.false

    def __init__(self, fmri_image, outputs=[], **keywords):
        LinearModelIterator.__init__(self, fmri_image, outputs, **keywords)
        
        self.fmri_image = fMRIImage(fmri_image)
        if self.normalize is not None:
            self.fmri_image.postread = self.normalize

        ftime = self.fmri_image.frametimes + self.tshift
        self.dmatrix = self.formula.design(time=ftime)

        if self.resid or self.output_fwhm:
            self.resid_output = ResidOutput(self.fmri_image.grid, path=self.path, basename='OLSresid', clobber=self.clobber)
            self.outputs.append(self.resid_output)

        self.rho_estimator = AROutput(self.fmri_image.grid, clobber=self.clobber)
        self.rho_estimator.setup_bias_correct(OLSModel(design=self.dmatrix))
        self.outputs.append(self.rho_estimator)

        self.setup_output()
        
    def model(self):
        ftime = self.fmri_image.frametimes + self.tshift
        if self.slicetimes is not None:
            _slice = self.iterator.grid.itervalue().slice
            model = OLSModel(design=self.formula.design(time=ftime + self.slicetimes[_slice[1]]))
        else:
            model = OLSModel(design=self.dmatrix)
        return model

    def fit(self, reference=None, **keywords):

        LinearModelIterator.fit(self, **keywords)

        sgrid = self.fmri_image.grid.subgrid(0)

        if self.output_fwhm:
            resid = fMRIImage(self.resid_output.img)
            fwhmest = fastFWHM(resid, fwhm=os.path.join(self.path, 'fwhmOLS.img'), clobber=self.clobber)
            fwhmest()
            self.fwhm_data = fwhmest.integrate(mask=self.mask)[1]
            print 'FWHM for data estimated as: %02f' % self.fwhm_data

        if reference is not None: self.estimateFWHM_AR(reference)
            
##      this will fail for non-affine grids, or grids
##      whose axes are not aligned in the standard way            

        #sigma = fwhm2sigma(self.fwhm_rho / N.array(self.fmri_image.image.pixdim[1:4][::-1]))

##         smoother = kernel_smooth.LinearFilter(sgrid, fwhm=self.fwhm_rho)
##         self.rho_estimator.img.grid = sgrid
##         self.rho = smoother.smooth(self.rho_estimator.img)

##        srho = scipy.ndimage.gaussian_filter(self.rho_estimator.img.readall(), sigma)
 ##       self.rho_estimator.img.writeslice(slice(0, self.rho_estimator.img.grid.shape[0], 1), srho)
        self.rho = self.rho_estimator.img
        self.getparcelmap()

    def getparcelmap(self):

        if self.mask is not None:
            _mask = self.mask.readall()
            self.rho[:] = N.where(_mask, self.rho[:], N.nan)

        if self.slicetimes == None:
            tmp = N.around(self.rho.readall() * (self.nmax / 2.)) / (self.nmax / 2.)
            tmp.shape = N.product(tmp.shape)
            parcelmap = tmp
            tmp = N.compress(1 - N.isnan(tmp), tmp)
            parcelseq = list(N.unique(tmp))
            self.fmri_image.grid.set_iter_param("parcelmap", parcelmap)
            self.fmri_image.grid.set_iter_param("parcelseq", parcelseq)            
            del(tmp); gc.collect()
        else:
            parcelmap = []
            parcelseq = []
            for i in range(self.rho.grid.shape[0]):
                tmp = self.rho[slice(i,i+1)]
                tmp.shape = N.product(tmp.shape)
                tmp = N.around(tmp * (self.nmax / 2.)) / (self.nmax / 2.)
                newlabels = list(N.unique(tmp))
                parcelseq += [newlabels]
                parcelmap += [tmp]
            self.fmri_image.grid.set_iter_param("parcelseq", parcelseq)
            self.fmri_image.grid.set_iter_param("parcelmap", parcelmap)            

    def setup_output(self):

        if not os.path.exists(self.path):
            os.makedirs(self.path)

        ftime = self.fmri_image.frametimes + self.tshift
        dmatrix = self.dmatrix.astype('<f8')
        ftime.shape = (ftime.shape[0],1)
        dmatrix = N.hstack([ftime, dmatrix])
        ftime.shape = (ftime.shape[0],)

        outname = os.path.join(self.path, 'matrix.csv')
        outfile = file(outname, 'w')
        tmatrix = [[fpformat.fix(dmatrix[i][j], 4) for j in range(dmatrix.shape[1])] for i in range(dmatrix.shape[0])]
        outfile.write('\n'.join(','.join(tmatrix[i]) for i in range(dmatrix.shape[0])))
        outfile.close()

        outname = os.path.join(self.path, 'matrix.bin')
        outfile = file(outname, 'w')
        dmatrix = self.dmatrix.astype('<f8')
        dmatrix.tofile(outfile)
        outfile.close()<commit_msg>cleanups on the RegrssionOutput class hierarchy, N.product(array.shape) -> array.size<commit_after>import gc, os, fpformat
from neuroimaging import traits

import numpy as N
import numpy.linalg as L
import numpy.random as R
import scipy.ndimage
from scipy.sandbox.models.utils import monotone_fn_inverter, rank 
from scipy.sandbox.models.regression import OLSModel, ARModel

from neuroimaging.modalities.fmri import fMRIImage
from neuroimaging.modalities.fmri.fmristat.delay import DelayContrast, DelayContrastOutput
from neuroimaging.modalities.fmri.regression import AROutput, TContrastOutput, \
  FContrastOutput, ResidOutput
from neuroimaging.algorithms import kernel_smooth
from neuroimaging.algorithms.fwhm import fastFWHM
from neuroimaging.algorithms.utils import fwhm2sigma
from neuroimaging.algorithms.statistics.regression import LinearModelIterator


# Imports related to pylab

from neuroimaging.defines import pylab_def
PYLAB_DEF, pylab = pylab_def()
if PYLAB_DEF:
    from neuroimaging.ui.visualization.multiplot import MultiPlot

class WholeBrainNormalize(traits.HasTraits):

    mask = traits.Any()

    def __init__(self, fmri_image, **keywords):
        traits.HasTraits.__init__(self, **keywords)
        if self.mask is not None:
            self._mask = self.mask.readall()
            self._mask.shape = self._mask.size
            
        self.n = fmri_image.grid.shape[0]
        self.avg = N.zeros((self.n,), N.float64)

        for i in range(self.n):
            d = fmri_image[slice(i,i+1)]
            if hasattr(self, '_mask'):
                d.shape = d.size
                d = N.compress(self._mask, d)
            self.avg[i] = d.mean()

    def __call__(self, fmri_data):
        out = N.zeros(fmri_data.shape, N.float64)
        for i in range(self.n):
            out[i] = fmri_data[i] * 100. / self.avg[i]
        return out

class fMRIStatOLS(LinearModelIterator):

    """
    OLS pass of fMRIstat.
    """
    
    normalize = traits.Any()
    
    formula = traits.Any()
    slicetimes = traits.Any()
    tshift = traits.Float(0.)

    fwhm_rho = traits.Float(6.)
    fwhm_data = traits.Float(6.)
    target_df = traits.Float(100.)
    nmax = traits.Int(200) # maximum number of rho values
    mask = traits.Any()
    path = traits.String('fmristat_run')
    resid = traits.false
    clobber = traits.false
    output_fwhm = traits.false

    def __init__(self, fmri_image, outputs=[], **keywords):
        LinearModelIterator.__init__(self, fmri_image, outputs, **keywords)
        
        self.fmri_image = fMRIImage(fmri_image)
        if self.normalize is not None:
            self.fmri_image.postread = self.normalize

        ftime = self.fmri_image.frametimes + self.tshift
        self.dmatrix = self.formula.design(time=ftime)

        if self.resid or self.output_fwhm:
            self.resid_output = ResidOutput(self.fmri_image.grid, path=self.path, basename='OLSresid', clobber=self.clobber)
            self.outputs.append(self.resid_output)

        model = OLSModel(design=self.dmatrix)
        self.rho_estimator = AROutput(self.fmri_image.grid, model, clobber=self.clobber)
        self.outputs.append(self.rho_estimator)

        self.setup_output()
        
    def model(self):
        ftime = self.fmri_image.frametimes + self.tshift
        if self.slicetimes is not None:
            _slice = self.iterator.grid.itervalue().slice
            model = OLSModel(design=self.formula.design(time=ftime + self.slicetimes[_slice[1]]))
        else:
            model = OLSModel(design=self.dmatrix)
        return model

    def fit(self, reference=None, **keywords):

        LinearModelIterator.fit(self, **keywords)

        sgrid = self.fmri_image.grid.subgrid(0)

        if self.output_fwhm:
            resid = fMRIImage(self.resid_output.img)
            fwhmest = fastFWHM(resid, fwhm=os.path.join(self.path, 'fwhmOLS.img'), clobber=self.clobber)
            fwhmest()
            self.fwhm_data = fwhmest.integrate(mask=self.mask)[1]
            print 'FWHM for data estimated as: %02f' % self.fwhm_data

        if reference is not None: self.estimateFWHM_AR(reference)
            
##      this will fail for non-affine grids, or grids
##      whose axes are not aligned in the standard way            

        #sigma = fwhm2sigma(self.fwhm_rho / N.array(self.fmri_image.image.pixdim[1:4][::-1]))

##         smoother = kernel_smooth.LinearFilter(sgrid, fwhm=self.fwhm_rho)
##         self.rho_estimator.img.grid = sgrid
##         self.rho = smoother.smooth(self.rho_estimator.img)

##        srho = scipy.ndimage.gaussian_filter(self.rho_estimator.img.readall(), sigma)
 ##       self.rho_estimator.img.writeslice(slice(0, self.rho_estimator.img.grid.shape[0], 1), srho)
        self.rho = self.rho_estimator.img
        self.getparcelmap()

    def getparcelmap(self):

        if self.mask is not None:
            _mask = self.mask.readall()
            self.rho[:] = N.where(_mask, self.rho[:], N.nan)

        if self.slicetimes == None:
            tmp = N.around(self.rho.readall() * (self.nmax / 2.)) / (self.nmax / 2.)
            tmp.shape = tmp.size
            parcelmap = tmp
            tmp = N.compress(1 - N.isnan(tmp), tmp)
            parcelseq = list(N.unique(tmp))
            self.fmri_image.grid.set_iter_param("parcelmap", parcelmap)
            self.fmri_image.grid.set_iter_param("parcelseq", parcelseq)            
            del(tmp); gc.collect()
        else:
            parcelmap = []
            parcelseq = []
            for i in range(self.rho.grid.shape[0]):
                tmp = self.rho[slice(i,i+1)]
                tmp.shape = tmp.size
                tmp = N.around(tmp * (self.nmax / 2.)) / (self.nmax / 2.)
                newlabels = list(N.unique(tmp))
                parcelseq += [newlabels]
                parcelmap += [tmp]
            self.fmri_image.grid.set_iter_param("parcelseq", parcelseq)
            self.fmri_image.grid.set_iter_param("parcelmap", parcelmap)            

    def setup_output(self):

        if not os.path.exists(self.path):
            os.makedirs(self.path)

        ftime = self.fmri_image.frametimes + self.tshift
        dmatrix = self.dmatrix.astype('<f8')
        ftime.shape = (ftime.shape[0],1)
        dmatrix = N.hstack([ftime, dmatrix])
        ftime.shape = (ftime.shape[0],)

        outname = os.path.join(self.path, 'matrix.csv')
        outfile = file(outname, 'w')
        tmatrix = [[fpformat.fix(dmatrix[i][j], 4) for j in range(dmatrix.shape[1])] for i in range(dmatrix.shape[0])]
        outfile.write('\n'.join(','.join(tmatrix[i]) for i in range(dmatrix.shape[0])))
        outfile.close()

        outname = os.path.join(self.path, 'matrix.bin')
        outfile = file(outname, 'w')
        dmatrix = self.dmatrix.astype('<f8')
        dmatrix.tofile(outfile)
        outfile.close()
{'fraction_non_alphanumeric': 0.11900946831755281, 'fraction_numerical': 0.006846321922796795, 'mean_word_length': 2.469176351692774, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6982607', 'n_tokens_mistral': 4813, 'n_tokens_neox': 4527, 'n_words': 848}