<commit_before>from nipype.interfaces.base import BaseInterface, TraitedSpec, traits,\
    InputMultiPath, File, OutputMultiPath
import numpy as np
import nibabel as nb
import nipy.neurospin.utils.design_matrix as dm
import nipy.neurospin.glm as GLM
import pylab
import os
from nipy.neurospin.mask import compute_mask
from nipype.utils.misc import isdefined

class FitGLMInputSpec(TraitedSpec):
    session_info = traits.Any(exists=True, desc='Session specific information generated by ``modelgen.SpecifyModel``')
    hrf_model = traits.Enum('Canonical', 'Canonical With Derivative', 'FIR', usedefault=True)
    drift_model = traits.Enum("Cosine", "Polynomial", "Blank", usedefault=True)
    functional_runs = InputMultiPath(File(exists=True), mandatory=True)
    TR = traits.Float(mandatory=True)
    model = traits.Enum("ar1", "spherical", usedefault=True)
    method = traits.Enum("kalman", "ols", usedefault=True)
    mask = traits.File(exists=True)
    
class FitGLMOutputSpec(TraitedSpec):
    beta = File(exists=True)
    nvbeta = traits.Any()
    s2 = File(exists=True)
    dof = traits.Any()
    constants = traits.Any()
    axis = traits.Any()
    reg_names = traits.List()
    residuals = traits.File(exists=True)
    
class FitGLM(BaseInterface):
    input_spec = FitGLMInputSpec
    output_spec = FitGLMOutputSpec
    
    def _run_interface(self, runtime):
        
        session_info = self.inputs.session_info
        
        nii = nb.load(self.inputs.functional_runs[0])              
        timeseries = nii.get_data()
        
        if isdefined(self.inputs.mask):
            mask = nb.load(self.inputs.mask).get_data() > 0
        else:
            mask = np.ones(timeseries.shape[:3]) == 1
            
        timeseries = timeseries[mask,:]
        
        for functional_run in self.inputs.functional_runs[1:]:
            nii = nb.load(functional_run)              
            timeseries = np.concatenate((timeseries,nii.get_data()[mask,:]), axis=1)
            
        nscans = timeseries.shape[1]
        
        if 'hpf' in session_info[0].keys():
            hpf = session_info[0]['hpf']
            drift_model=self.inputs.drift_model
        else:
            hpf=0
            drift_model = "Blank"
        
        reg_names = []
        for reg in session_info[0]['regress']:
            reg_names.append(reg['name'])
        
        reg_vals = np.zeros((nscans,len(reg_names)))
        for i in range(len(reg_names)):
            reg_vals[:,i] = np.array(session_info[0]['regress'][i]['val']).reshape(1,-1)
        
        
        frametimes= np.linspace(0, (nscans-1)*self.inputs.TR, nscans)
        
        conditions = []
        onsets = []
        duration = []
        
        for i,cond in enumerate(session_info[0]['cond']):
            onsets += cond['onset']
            conditions += [cond['name']]*len(cond['onset'])
            if len(cond['duration']) == 1:
                duration += cond['duration']*len(cond['onset'])
                
        
        paradigm =  dm.BlockParadigm(con_id=conditions, onset=onsets, duration=duration)
        design_matrix, self._reg_names = dm.dmtx_light(frametimes, paradigm, drift_model=drift_model, hfcut=hpf,
               hrf_model=self.inputs.hrf_model, 
               add_regs=reg_vals,
               add_reg_names=reg_names
               )
        pylab.pcolor(design_matrix)
        pylab.savefig("design_matrix.pdf")
        
        glm = GLM.glm()
        glm.fit(timeseries.T, design_matrix, method=self.inputs.method, model=self.inputs.model)
        
        self._beta_file = os.path.abspath("beta.nii")
        beta = np.zeros(mask.shape + (glm.beta.shape[0],))
        beta[mask,:] = glm.beta.T
        nb.save(nb.Nifti1Image(beta, nii.get_affine()), self._beta_file)
        
        self._s2_file = os.path.abspath("s2.nii")
        s2 = np.zeros(mask.shape)
        s2[mask] = glm.s2
        nb.save(nb.Nifti1Image(s2, nii.get_affine()), self._s2_file)
        <commit_msg>Added design matrix normalization<commit_after>from nipype.interfaces.base import BaseInterface, TraitedSpec, traits,\
    InputMultiPath, File, OutputMultiPath
import numpy as np
import nibabel as nb
import nipy.neurospin.utils.design_matrix as dm
import nipy.neurospin.glm as GLM
import pylab
import os
from nipy.neurospin.mask import compute_mask
from nipype.utils.misc import isdefined

class FitGLMInputSpec(TraitedSpec):
    session_info = traits.Any(exists=True, desc='Session specific information generated by ``modelgen.SpecifyModel``')
    hrf_model = traits.Enum('Canonical', 'Canonical With Derivative', 'FIR', usedefault=True)
    drift_model = traits.Enum("Cosine", "Polynomial", "Blank", usedefault=True)
    functional_runs = InputMultiPath(File(exists=True), mandatory=True)
    TR = traits.Float(mandatory=True)
    model = traits.Enum("ar1", "spherical", usedefault=True)
    method = traits.Enum("kalman", "ols", usedefault=True)
    mask = traits.File(exists=True)
    normalize_design_matrix = traits.Bool(False, usedefault=True)
    
class FitGLMOutputSpec(TraitedSpec):
    beta = File(exists=True)
    nvbeta = traits.Any()
    s2 = File(exists=True)
    dof = traits.Any()
    constants = traits.Any()
    axis = traits.Any()
    reg_names = traits.List()
    residuals = traits.File(exists=True)
    
class FitGLM(BaseInterface):
    input_spec = FitGLMInputSpec
    output_spec = FitGLMOutputSpec
    
    def _run_interface(self, runtime):
        
        session_info = self.inputs.session_info
        
        nii = nb.load(self.inputs.functional_runs[0])              
        timeseries = nii.get_data()
        
        if isdefined(self.inputs.mask):
            mask = nb.load(self.inputs.mask).get_data() > 0
        else:
            mask = np.ones(timeseries.shape[:3]) == 1
            
        timeseries = timeseries[mask,:]
        
        for functional_run in self.inputs.functional_runs[1:]:
            nii = nb.load(functional_run)              
            timeseries = np.concatenate((timeseries,nii.get_data()[mask,:]), axis=1)
            
        nscans = timeseries.shape[1]
        
        if 'hpf' in session_info[0].keys():
            hpf = session_info[0]['hpf']
            drift_model=self.inputs.drift_model
        else:
            hpf=0
            drift_model = "Blank"
        
        reg_names = []
        for reg in session_info[0]['regress']:
            reg_names.append(reg['name'])
        
        reg_vals = np.zeros((nscans,len(reg_names)))
        for i in range(len(reg_names)):
            reg_vals[:,i] = np.array(session_info[0]['regress'][i]['val']).reshape(1,-1)
        
        
        frametimes= np.linspace(0, (nscans-1)*self.inputs.TR, nscans)
        
        conditions = []
        onsets = []
        duration = []
        
        for i,cond in enumerate(session_info[0]['cond']):
            onsets += cond['onset']
            conditions += [cond['name']]*len(cond['onset'])
            if len(cond['duration']) == 1:
                duration += cond['duration']*len(cond['onset'])
                
        
        paradigm =  dm.BlockParadigm(con_id=conditions, onset=onsets, duration=duration)
        design_matrix, self._reg_names = dm.dmtx_light(frametimes, paradigm, drift_model=drift_model, hfcut=hpf,
               hrf_model=self.inputs.hrf_model, 
               add_regs=reg_vals,
               add_reg_names=reg_names
               )
        print design_matrix.shape
        print len(self._reg_names)
        if self.inputs.normalize_design_matrix:
            for i in range(len(self._reg_names)-1):
                design_matrix[:,i] = (design_matrix[:,i]-design_matrix[:,i].mean())/design_matrix[:,i].std()
        pylab.pcolor(design_matrix)
        pylab.savefig("design_matrix.pdf")
        
        glm = GLM.glm()
        glm.fit(timeseries.T, design_matrix, method=self.inputs.method, model=self.inputs.model)
        
        self._beta_file = os.path.abspath("beta.nii")
        beta = np.zeros(mask.shape + (glm.beta.shape[0],))
        beta[mask,:] = glm.beta.T
        nb.save(nb.Nifti1Image(beta, nii.get_affine()), self._beta_file)
        
        self._s2_file = os.path.abspath("s2.nii")
        s2 = np.zeros(mask.shape)
        s2[mask] = glm.s2
        nb.save(nb.Nifti1Image(s2, nii.get_affine()), self._s2_file)
        
{'fraction_non_alphanumeric': 0.12613582018173122, 'fraction_numerical': 0.007293161166905786, 'mean_word_length': 2.42827868852459, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, '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': '2623493', 'n_tokens_mistral': 2791, 'n_tokens_neox': 2739, 'n_words': 446}