Python pyhrf.verbose函数代码示例

    def test_gls_recursive(self):
        cmd = ["pyhrf_gls", "-r", self.tmp_dir]
        output = check_output(cmd)

        pyhrf.verbose(1, "output:")
        pyhrf.verbose(1, output)
        expected_ouput = """%s:
%s/subject1:
%s/subject1/fmri:
paradigm.csv

%s/subject1/fmri/analysis:
analysis_result_1.nii
analysis_result_2.csv
analysis_summary.txt

%s/subject1/fmri/run1:
bold_scan_[1...3].nii

%s/subject1/fmri/run2:
bold_scan_[1...3].nii

%s/subject1/t1mri:
anatomy.{hdr,img}

""" % (
            (self.tmp_dir,) * 7
        )
        if output != expected_ouput:
            raise Exception(
                "Output of command %s is not as expected.\n"
                "Output is:\n%sExcepted:\n%s" % (" ".join(cmd), output, expected_ouput)
            )

    def __init__(self, sampled_variables, nb_its_max, obs_pace=1, burnin=.3,
                 sample_hist_pace=-1, obs_hist_pace=-1,):

        self.variables = {}
        self.sampled_variables = sampled_variables

        for v in sampled_variables:
            self.set_variable(v.name, v)

        def get_fraction_or_nb(nb, tot):
            if nb>0. and nb<1.:
                return int(round(tot * nb))
            else:
                return nb

        self.nb_its_max = nb_its_max
        self.burnin = get_fraction_or_nb(burnin, nb_its_max)
        self.smpl_hist_pace = get_fraction_or_nb(sample_hist_pace, nb_its_max)
        self.obs_hist_pace = get_fraction_or_nb(obs_hist_pace, nb_its_max)
        self.tracked_quantities = {}

        pyhrf.verbose(1, 'GibbsSampler init. Burnin: %d, nb_its_max: %d, '\
                      'smpl_hist_pace: %d, obs_hist_pace: %d,'\
                      %(self.burnin, self.nb_its_max, self.smpl_hist_pace,
                        self.obs_hist_pace))

    def computeComponentsApost(self, variables, j, XhtQXh):
        sIMixtP = variables[self.samplerEngine.I_MIXT_PARAM]
        var = sIMixtP.getCurrentVars()
        mean = sIMixtP.getCurrentMeans()
        rb = variables[self.samplerEngine.I_NOISE_VAR].currentValue
        nrls = self.currentValue
        
        #for j in xrange(self.nbConditions):
        gTQgjrb = XhtQXh[:,j]/rb   # de taille nbVox
    
        ej = self.varYtilde + repmat(nrls[j,:], self.ny, 1) * self.varXh[:,:,j].swapaxes(0,1)
        numpy.divide(diag(dot(self.varXhtQ[:,j,:],ej)), rb, self.varXjhtQjeji)

        for c in xrange(self.nbClasses):            #ici classe: 0 (inactif) ou 1 (actif)
            self.varClassApost[c,j,:] = 1./(1./var[c,j] + gTQgjrb)
            numpy.sqrt(self.varClassApost[c,j,:], self.sigClassApost[c,j,:])
            if c > 0: # assume 0 stands for inactivating class
                numpy.multiply(self.varClassApost[c,j,:],
                               add(mean[c,j]/var[c,j], self.varXjhtQjeji),
                               self.meanClassApost[c,j,:])
            else:
                multiply(self.varClassApost[c,j,:], self.varXjhtQjeji,
                         self.meanClassApost[c,j,:])
                
            pyhrf.verbose(5, 'meanClassApost %d cond %d :'%(c,j))
            pyhrf.verbose.printNdarray(5, self.meanClassApost[c,j,:])

    def packSamplerInput(self, roiData):

        try:
            shrf = self.sampler.getVariable('hrf')
        except KeyError:
            shrf = self.sampler.getVariable('brf')
            
        hrfDuration = shrf.duration
        zc = shrf.zc

        simu = None

        if simu != None and shrf.sampleFlag==0:
            hrfDuration = (len(simu.hrf.get_hrf(0,0))-1)*simu.hrf.dt
            pyhrf.verbose(6,'Found simulation data and hrf is '\
                          'not sampled, setting hrfDuration to:' \
                          +str(hrfDuration))

        pyhrf.verbose(2,'building BOLDSamplerInput ...')

        if simu == None or shrf.sampleFlag:
            dt = self.dt if (self.dt!=None and self.dt!=0.) else -self.dtMin
        elif simu != None and shrf.sampleFlag == 0:
            dt = simu.hrf.dt

        samplerInput = self.sampler.inputClass(roiData, dt=dt,
                                               typeLFD=self.driftLfdType,
                                               paramLFD=self.driftLfdParam,
                                               hrfZc=zc,
                                               hrfDuration=hrfDuration)
        return samplerInput

 def _compute_graph(self):
     if self.data_type != 'volume':
         raise Exception('Can only compute graph for volume data')
     pyhrf.verbose(6, 'FmriData._compute_graph() ...')
     to_discard = [self.backgroundLabel]
     self._graph = parcels_to_graphs(self.roiMask, kerMask3D_6n,
                                     toDiscard=to_discard)

    def checkAndSetInitValue(self, variables):
        smplVarDrift = variables[self.samplerEngine.I_ETA]
        smplVarDrift.checkAndSetInitValue(variables)
        varDrift = smplVarDrift.currentValue

        if self.useTrueValue :
            if self.trueValue is not None:
                self.currentValue = self.trueValue
            else:
                raise Exception('Needed a true value for drift init but '\
                                    'None defined')

        if 0 and self.currentValue is None :
            #if not self.sampleFlag and self.dataInput.simulData is None :
                #self.currentValue = self.dataInput.simulData.drift.lfd
                #pyhrf.verbose(6,'drift dimensions :' \
                              #+str(self.currentValue.shape))
                #pyhrf.verbose(6,'self.dimDrift :' \
                              #+str(self.dimDrift))
                #assert self.dimDrift == self.currentValue.shape[0]
            #else:
            self.currentValue = np.sqrt(varDrift) * \
                np.random.randn(self.dimDrift, self.nbVox)

        if self.currentValue is None:
            pyhrf.verbose(1,"Initialisation of Drift from the data")
            ptp = numpy.dot(self.P.transpose(),self.P)
            invptp = numpy.linalg.inv(ptp)
            invptppt = numpy.dot(invptp, self.P.transpose())
            self.currentValue = numpy.dot(invptppt,self.dataInput.varMBY)

        self.updateNorm()
        self.matPl = dot(self.P, self.currentValue)
        self.ones_Q_J = np.ones((self.dimDrift, self.nbVox))
        self.ones_Q   = np.ones((self.dimDrift))

def remote_copy(files, host, user, path, transfer_tool='ssh'):
    if transfer_tool == 'paramiko':
        import paramiko
        pyhrf.verbose(1, 'Copying files to remote destination %[email protected]%s:%s ...' \
                          %(host,user,path))
        ssh = paramiko.SSHClient()
        known_hosts_file = os.path.join("~", ".ssh", "known_hosts")
        ssh.load_host_keys(os.path.expanduser(known_hosts_file))
        ssh.connect(host, username=user)
        sftp = ssh.open_sftp()
        for f in files:
            remotepath = op.join(path,op.basename(f))
            pyhrf.verbose(2, f + ' -> ' + remotepath + ' ...')
            flocal = open(f)
            remote_file = sftp.file(remotepath, "wb")
            remote_file.set_pipelined(True)
            remote_file.write(flocal.read())
            flocal.close()
            remote_file.close()
        sftp.close()
        ssh.close()
    else:
        sfiles = string.join(['"%s"'%f for f in files], ' ')

        scp_cmd = 'scp -C %s "%[email protected]%s:%s"' %(sfiles, user, host, path)
        pyhrf.verbose(1, 'Data files transfer with scp ...')
        pyhrf.verbose(2, scp_cmd)
        if os.system(scp_cmd) != 0:
            raise Exception('Error while scp ...')

    pyhrf.verbose(1, 'Copy done!')

    return [op.join(path,op.basename(f)) for f in files]

    def sampleNextInternal(self, variables):
        #TODO : comment
        sIMixtP = variables[self.samplerEngine.I_MIXT_PARAM_NRLS_BAR]
        varCI = sIMixtP.currentValue[sIMixtP.I_VAR_CI]
        varCA = sIMixtP.currentValue[sIMixtP.I_VAR_CA]
        meanCA = sIMixtP.currentValue[sIMixtP.I_MEAN_CA]

        self.labelsSamples = rand(self.nbConditions, self.nbVox)
        self.nrlsSamples = randn(self.nbConditions, self.nbVox)

        if self.imm:
            #self.sampleNrlsParallel(varXh, rb, h, varLambda, varCI,
            #                        varCA, meanCA, gTQg, variables)
            raise NotImplementedError("IMM with drift sampling is not available")
        else: #smm
            self.sampleNrlsSerial(varCI, varCA, meanCA, variables)
            #self.computeVarYTildeOpt(varXh)
            #matPl = self.samplerEngine.getVariable('drift').matPl
            #self.varYbar = self.varYtilde - matPl

        if (self.currentValue >= 1000).any() and \
                pyhrf.__usemode__ == pyhrf.DEVEL:
            pyhrf.verbose(2, "Weird NRL values detected ! %d/%d" \
                              %((self.currentValue >= 1000).sum(),
                                self.nbVox*self.nbConditions) )
            #pyhrf.verbose.set_verbosity(6)

        if pyhrf.verbose.verbosity >= 4:
            self.reportDetection()

        self.printState(4)
        self.iteration += 1 #TODO : factorize !!

    def from_py_object(self, label, obj, parent=None):

        pyhrf.verbose(6, "UiNode.from_py_object(label=%s,obj=%s) ..." % (label, str(obj)))

        if isinstance(obj, Initable):
            n = obj.to_ui_node(label, parent)
        else:
            if UiNode._pyobj_has_leaf_type(obj):
                if isinstance(obj, np.ndarray):
                    dt = str(obj.dtype.name)
                    sh = str(obj.shape)
                    n = UiNode(label, attributes={"type": "ndarray", "dtype": dt, "shape": sh})
                    s = " ".join(str(e) for e in obj.ravel())
                    n.add_child(UiNode(s))
                elif obj is None:
                    n = UiNode(label, attributes={"type": "None"})
                    n.add_child(UiNode("None"))
                else:
                    n = UiNode(label, attributes={"type": obj.__class__})
                    n.add_child(UiNode(str(obj)))
            elif isinstance(obj, list):
                n = UiNode(label, attributes={"type": "list"})
                for i, sub_val in enumerate(obj):
                    n.add_child(UiNode.from_py_object("item%d" % i, sub_val))
            elif isinstance(obj, (dict, OrderedDict)):
                t = ["odict", "dict"][obj.__class__ == dict]
                n = UiNode(label, attributes={"type": t})
                for k, v in obj.iteritems():
                    n.add_child(UiNode.from_py_object(k, v))
            else:
                raise Exception(
                    "In UiNode.from_py_object, unsupported object: " "%s (type: %s)" % (str(obj), str(type(obj)))
                )
        return n

    def analyse_roi(self, atomData):
        """
        Launch the JDE Gibbs Sampler on a parcel-specific data set *atomData*
        Args:
            - atomData (pyhrf.core.FmriData): parcel-specific data
        Returns:
            JDE sampler object
        """
        #print 'atomData:', atomData

        if self.copy_sampler:
            sampler = copyModule.deepcopy(self.sampler)
        else:
            sampler = self.sampler
        sInput = self.packSamplerInput(atomData)
        sampler.linkToData(sInput)
        #if self.parameters[self.P_RANDOM_SEED] is not None:
        #    np.random.seed(self.parameters[self.P_RANDOM_SEED])
        # #HACK:
        # if len(self.roi_ids) > 0:
        #     if atomData.get_roi_id() not in self.roi_ids:
        #         return None

        pyhrf.verbose(1, 'Treating region %d' %(atomData.get_roi_id()))
        sampler.runSampling()
        pyhrf.verbose(1, 'Cleaning memory ...')
        sampler.dataInput.cleanMem()
        return sampler

def project_fmri(input_mesh, data_file, output_tex_file, 
                 output_kernels_file=None, data_resolution=None,
                 geod_decay=5., norm_decay=2., kernel_size=7, 
                 tex_bin_threshold=None):

    if output_kernels_file is None:
        tmp_dir = tempfile.mkdtemp(prefix='pyhrf_surf_proj',
                                   dir=pyhrf.cfg['global']['tmp_path'])

        kernels_file = op.join(tmp_dir, add_suffix(op.basename(data_file),
                                                   '_kernels'))
        tmp_kernels_file = True
    else:
        kernels_file = output_kernels_file
        tmp_kernels_file = False

    if data_resolution is not None:
        resolution = data_resolution
    else:
        resolution = read_spatial_resolution(data_file)

    pyhrf.verbose(1,'Data resolution: %s' %resolution)
    pyhrf.verbose(2,'Projection parameters:')
    pyhrf.verbose(2,'   - geodesic decay: %f mm' %geod_decay)
    pyhrf.verbose(2,'   - normal decay: %f mm' %norm_decay)
    pyhrf.verbose(2,'   - kernel size: %f voxels' %kernel_size)
    
    create_projection_kernels(input_mesh, kernels_file, resolution, 
                              geod_decay, norm_decay, kernel_size)
    
    project_fmri_from_kernels(input_mesh, kernels_file, data_file,
                              output_tex_file, tex_bin_threshold)

    if tmp_kernels_file:
        os.remove(kernels_file)

    def test_ward_spatial_real_data(self):
        from pyhrf.glm import glm_nipy_from_files
        #pyhrf.verbose.verbosity = 2

        fn = 'subj0_parcellation.nii.gz'
        mask_file = pyhrf.get_data_file_name(fn)

        bold = 'subj0_bold_session0.nii.gz'
        bold_file = pyhrf.get_data_file_name(bold)

        paradigm_csv_file = pyhrf.get_data_file_name('paradigm_loc_av.csv')
        output_dir = self.tmp_dir
        output_file = op.join(output_dir,
                              'parcellation_output_test_real_data.nii')

        tr=2.4
        bet = glm_nipy_from_files(bold_file, tr,
                                paradigm_csv_file, output_dir,
                                mask_file, session=0,
                                contrasts=None,
                                hrf_model='Canonical',
                                drift_model='Cosine', hfcut=128,
                                residuals_model='spherical',
                                fit_method='ols', fir_delays=[0])[0]

        pyhrf.verbose(2, 'betas_files: %s' %' '.join(bet))

        cmd = 'pyhrf_parcellate_glm -m %s %s -o %s -v %d -n %d '\
            '-t ward_spatial ' \
        %(mask_file, ' '.join(bet), output_file,
          pyhrf.verbose.verbosity, 10)

        if os.system(cmd) != 0 :
            raise Exception('"' + cmd + '" did not execute correctly')
        pyhrf.verbose(1, 'cmd: %s' %cmd)

    def checkAndSetInitHabit(self, variables) :
     
        # init habituation speed factors and time-varying NRLs
        if self.habits == None :  # if no habituation speed specified 
            if not self.sampleHabitFlag :
                if self.dataInput.simulData != None : # Attention: on a un probleme si on fait tourner ce modele sur des donnees simulees par le modele stationnaire. Dans ce cas, il faut forcer ici a passer et prendre des habits nulles
                   ## using simulated Data for HABITUATION sampling
                   print "load Habituation from simulData", self.dataInput.simulData.habitspeed.data
                   self.habits = self.dataInput.simulData.habitspeed.data
                   self.timeNrls = self.dataInput.simulData.nrls.timeNrls
                   self.Gamma = self.setupGamma()
                else : # sinon, on prend que des zeros (modele stationnaire)
                    self.habits = numpy.zeros((self.nbConditions, self.nbVox), dtype=float)
                    self.setupTimeNrls()
            else:
                pyhrf.verbose(2, "Uniform set up of habituation factors")
                habitCurrent = numpy.zeros((self.nbConditions, self.nbVox), dtype=float)
                for nc in xrange(self.nbConditions):
                    habitCurrent[nc,self.voxIdx[1][nc]] = numpy.random.rand(self.cardClass[1][nc])
                self.habits = habitCurrent

        if self.outputRatio :
            self.ratio = zeros((self.nbConditions, self.nbVox, 2), dtype = float)
            self.ratiocourbe = zeros((self.nbConditions, self.nbVox, 100, 5), dtype = float)
            self.compteur = numpy.zeros((self.nbConditions, self.nbVox), dtype=float)
         
        self.setupTimeNrls()
        pyhrf.verbose(4, 'habituation initiale')
        pyhrf.verbose.printNdarray(5, self.habits)

def create_hrf_from_territories(hrf_territories, primary_hrfs):

    pyhrf.verbose(1,'create_hrf_from_territories ...')
    pyhrf.verbose(1,' inputs: hrf_territories %s, primary_hrfs (%d,%d)' \
                      %(str(hrf_territories.shape), len(primary_hrfs),
                        len(primary_hrfs[0][0])))
    assert hrf_territories.ndim == 1
    hrfs = np.zeros((hrf_territories.size,primary_hrfs[0][0].size))
    territories = np.unique(hrf_territories)
    territories.sort()
    if territories.min() == 1:
        territories = territories - 1

    assert territories.min() == 0
    assert territories.max() <= len(primary_hrfs)
    #print hrfs.shape

    #sm = ','.join(['m[%d]'%d for d in range(hrf_territories.ndim)] + [':'])
    for territory in territories:
        #TODO: test consitency in hrf lengths
        m = np.where(hrf_territories==territory)[0]
        # print 'm:',m
        # print 'hrfs[m,:].shape:', hrfs[m,:].shape
        # print 'primary_hrfs[territory][1].shape:', \
        #     primary_hrfs[territory][1].shape
        # print primary_hrfs[territory][1]
        #print hrfs[m,:].shape
        #exec('hrfs[%s] = primary_hrfs[territory][1]' %sm)
        hrfs[m,:] = primary_hrfs[territory][1]
        #print hrfs[m,:]
    return hrfs.transpose()

def create_asl_from_stim_induced(bold_stim_induced, perf_stim_induced,
                                 ctrl_tag_mat, dsf,
                                 perf_baseline, noise, drift=None, outliers=None):
    """
    Downsample stim_induced signal according to downsampling factor 'dsf' and
    add noise and drift (nuisance signals) which has to be at downsampled
    temporal resolution.
    """
    bold = bold_stim_induced[0:-1:dsf,:].copy()
    perf = np.dot(ctrl_tag_mat, (perf_stim_induced[0:-1:dsf,:].copy() + \
                                 perf_baseline))

    pyhrf.verbose(3, 'create_asl_from_stim_induced ...')
    pyhrf.verbose(3, 'bold shape: %s, perf_shape: %s, noise shape: %s, '\
                  'drift shape: %s' %(str(bold.shape), str(perf.shape),
                                      str(noise.shape), str(drift.shape)))

    asl = bold + perf
    if drift is not None:
        asl += drift
    if outliers is not None:
        asl += outliers
    asl += noise

    return asl