<commit_before>"""
The main routine of this package that aims at performing the
extraction of ROIs from multisubject dataset using the localization.
This has been puclished in Thirion et al. Structural Analysis of fMRI
Data Revisited: Improving the Sensitivity and Reliability of fMRI
Group Studies.  IEEE TMI 2007

Author : <NAME>, 2006-2008
"""

#autoindent 

import numpy as np
import numpy.random as nr


# import nipy.neurospin.clustering.clustering as fc
import nipy.neurospin.graph.BPmatch as BPmatch
import nipy.neurospin.graph.field as ff
import nipy.neurospin.graph.graph as fg
from nipy.neurospin.spatial_models import hroi 
from nipy.neurospin.clustering.hierarchical_clustering import \
    Average_Link_Graph_segment

class landmark_regions(hroi.NROI):
    """
    This class is intended to represent a set of inter-subject regions
    besides the standard multiple ROI features, its has
    features 'subjects_ids' and 'individual_positions' to describe
    from which subjects at which position it is found.
    """
    def __init__(self,k,parents=None,header=None,id=None,subj=None,coord=None):
        """
        Building the landmark_region

        INPUT
        - k (int): number of nodes/landmarks included into it
        nb :
        - parents = None: array of shape(self.k) describing the
        hierarchical relationship
        if None parents = np.arange(k) is sued instead
        - header (temporary): space-defining image header
        to embed the structure in an image space
        - subj=None: k-length list of subjects
        (these correspond to ROI feature)
        - coord=None; k-length list of  coordinate arrays
        CAVEAT:
        discrete is used here for subject identification
        """
        k = int(k)
        if k<1: raise ValueError, "cannot create an empty LR"
        if parents==None:
            parents = np.arange(k)
        hroi.NROI.__init__(self,parents,header,discrete=subj,id=id)
        self.set_discrete_feature('position',coord)
        self.subj = self.discrete

    def centers(self):
        """
        c = self.centers()
        returns the average of the coordinates for each region
        """
        centers = self.discrete_to_roi_features('position')
        return centers

    def homogeneity(self):
        """ returns the mean distance between points within each LR
        """
        import nipy.neurospin.eda.dimension_reduction  as dr 
        coord = self.discrete_features['position']
        size = self.get_size()
        h = np.zeros(self.k)
        for k in range(self.k):
             edk = dr.Euclidian_distance(coord[k]) 
             h[k] = edk.sum()/(size[k]*(size[k]-1))
        return h
             
    def HPD(self,k,cs,pval = 0.95,dmax=1.0):
        """
        hpd = self.HPD(k,cs,dmax = 10,pval = 0.95)
        Sample the postreior density of being in k
        on a grid defined by cs, assuming that the roi is an ellipsoid
        INPUT:
        - cs: an array of shape(n,dim) a set of input coordinates
         - pval=0.95 cutoff for the CR
         - dmax=1.0 : an upper bound for the spatial variance
        to avoid degenerate variance
        OUPUT:
        - hpd array of shape(n) that yields the value
        """
        if k>self.k:
            raise ValueError, 'wrong region index'
        
        coord = self.discrete_features['position'][k]
        centers = self.discrete_to_roi_features('position')
        dim = centers.shape[1]
        
        if cs.shape[1]!=dim:
            raise ValueError, "incompatible dimensions"
        
        dx = coord-centers[k]
        covariance = np.dot(dx.T,dx)/coord.shape[0]
        import numpy.linalg as L
        U,S,V = L.svd(covariance,0)
        eps = (dmax**2)/coord.shape[0]
        sqrts = np.sqrt(1/np.maximum(S,eps))
        dx = cs-centers[k]
        dx = np.dot(dx,U)
        dx = np.dot(dx,np.diag(sqrts))
        delta = np.sum(dx**2,1)
        lcst = -np.log(2*np.pi)*dim/2+(np.log(sqrts)).sum()
        hpd = np.exp(lcst-delta/2)

        import scipy.special as sp
        gamma = 2*sp.erfinv(pval)**2
        #
        #--- all the following is to solve the equation
        #--- erf(x/sqrt(2))-x*exp(-x**2/2)/sqrt(pi/2) = alpha
        #--- should better be put elsewhere
        #
        def dicho_solve_lfunc(alpha,eps = 1.e-7):
            if alpha>1:
                raise ValueError, "no solution for alpha>1"
            if alpha>1-1.e-15:
                return np.infty
            if alpha<0:
                raise ValueError, "no solution for alpha<0" 
            if alpha<1.e-15:
                return 0
    
            xmin = sp.erfinv(alpha)*np.sqrt(2)
            xmax = 2*xmin
            while lfunc(xmax)<alpha:
                xmax*=2
                xmin*=2
            return (dichomain_lfunc(xmin,xmax,eps,alpha))

        def dichomain_lfunc(xmin,xmax,eps,alpha):
            x =  (xmin+xmax)/2
            if xmax<xmin+eps:
                return x
            else:
                print xmin,x,xmax
                if lfunc(x)>alpha:
                    return dichomain_lfunc(xmin,x,eps,alpha)
                else:
                    return dichomain_lfunc(x,xmax,eps,alpha)
        
        def lfunc(x):
            return sp.erf(x/np.sqrt(2))-x*np.exp(-x**2/2)/np.sqrt(np.pi/2)

        gamma = dicho_solve_lfunc(pval)**2
        hpd[delta>gamma]=0
         return hpd
<commit_msg>a slight improvement of region definition<commit_after>"""
The main routine of this package that aims at performing the
extraction of ROIs from multisubject dataset using the localization.
This has been puclished in Thirion et al. Structural Analysis of fMRI
Data Revisited: Improving the Sensitivity and Reliability of fMRI
Group Studies.  IEEE TMI 2007

Author : <NAME>, 2006-2008
"""

#autoindent 

import numpy as np
import numpy.random as nr


# import nipy.neurospin.clustering.clustering as fc
import nipy.neurospin.graph.BPmatch as BPmatch
import nipy.neurospin.graph.field as ff
import nipy.neurospin.graph.graph as fg
from nipy.neurospin.spatial_models import hroi 
from nipy.neurospin.clustering.hierarchical_clustering import \
    Average_Link_Graph_segment

class landmark_regions(hroi.NROI):
    """
    This class is intended to represent a set of inter-subject regions
    besides the standard multiple ROI features, its has
    features 'subjects_ids' and 'individual_positions' to describe
    from which subjects at which position it is found.
    """
    def __init__(self,k,parents=None,header=None,id=None,subj=None,coord=None):
        """
        Building the landmark_region

        INPUT
        - k (int): number of nodes/landmarks included into it
        nb :
        - parents = None: array of shape(self.k) describing the
        hierarchical relationship
        if None parents = np.arange(k) is sued instead
        - header (temporary): space-defining image header
        to embed the structure in an image space
        - subj=None: k-length list of subjects
        (these correspond to ROI feature)
        - coord=None; k-length list of  coordinate arrays
        CAVEAT:
        discrete is used here for subject identification
        """
        k = int(k)
        if k<1: raise ValueError, "cannot create an empty LR"
        if parents==None:
            parents = np.arange(k)
        hroi.NROI.__init__(self,parents,header,discrete=subj,id=id)
        self.set_discrete_feature('position',coord)
        self.subj = self.discrete

    def centers(self):
        """
        c = self.centers()
        returns the average of the coordinates for each region
        """
        centers = self.discrete_to_roi_features('position')
        return centers

    def homogeneity(self):
        """ returns the mean distance between points within each LR
        """
        import nipy.neurospin.eda.dimension_reduction  as dr 
        coord = self.discrete_features['position']
        size = self.get_size()
        h = np.zeros(self.k)
        for k in range(self.k):
             edk = dr.Euclidian_distance(coord[k]) 
             h[k] = edk.sum()/(size[k]*(size[k]-1))
        return h
             
    def HPD(self,k,cs,pval = 0.95,dmax=1.0):
        """
        hpd = self.HPD(k,cs,dmax = 10,pval = 0.95)
        Sample the postreior density of being in k
        on a grid defined by cs, assuming that the roi is an ellipsoid
        INPUT:
        - cs: an array of shape(n,dim) a set of input coordinates
         - pval=0.95 cutoff for the CR
         - dmax=1.0 : an upper bound for the spatial variance
        to avoid degenerate variance
        OUPUT:
        - hpd array of shape(n) that yields the value
        """
        if k>self.k:
            raise ValueError, 'wrong region index'
        
        coord = self.discrete_features['position'][k]
        centers = self.discrete_to_roi_features('position')
        dim = centers.shape[1]
        
        if cs.shape[1]!=dim:
            raise ValueError, "incompatible dimensions"
        
        dx = coord-centers[k]
        covariance = np.dot(dx.T,dx)/coord.shape[0]
        import numpy.linalg as L
        U,S,V = L.svd(covariance,0)
        eps = (dmax**2)/coord.shape[0]
        sqrts = np.sqrt(1/np.maximum(S,eps))
        dx = cs-centers[k]
        dx = np.dot(dx,U)
        dx = np.dot(dx,np.diag(sqrts))
        delta = np.sum(dx**2,1)
        lcst = -np.log(2*np.pi)*dim/2+(np.log(sqrts)).sum()
        hpd = np.exp(lcst-delta/2)

        import scipy.special as sp
        gamma = 2*sp.erfinv(pval)**2
        #
        #--- all the following is to solve the equation
        #--- erf(x/sqrt(2))-x*exp(-x**2/2)/sqrt(pi/2) = alpha
        #--- should better be put elsewhere
        #
        def dicho_solve_lfunc(alpha,eps = 1.e-7):
            if alpha>1:
                raise ValueError, "no solution for alpha>1"
            if alpha>1-1.e-15:
                return np.infty
            if alpha<0:
                raise ValueError, "no solution for alpha<0" 
            if alpha<1.e-15:
                return 0
    
            xmin = sp.erfinv(alpha)*np.sqrt(2)
            xmax = 2*xmin
            while lfunc(xmax)<alpha:
                xmax*=2
                xmin*=2
            return (dichomain_lfunc(xmin,xmax,eps,alpha))

        def dichomain_lfunc(xmin,xmax,eps,alpha):
            x =  (xmin+xmax)/2
            if xmax<xmin+eps:
                return x
            else:
                print xmin,x,xmax
                if lfunc(x)>alpha:
                    return dichomain_lfunc(xmin,x,eps,alpha)
                else:
                    return dichomain_lfunc(x,xmax,eps,alpha)
        
        def lfunc(x):
            return sp.erf(x/np.sqrt(2))-x*np.exp(-x**2/2)/np.sqrt(np.pi/2)

        gamma = dicho_solve_lfunc(pval)**2
        hpd[delta>gamma]=0
        return hpd

{'fraction_non_alphanumeric': 0.09975108324882456, 'fraction_numerical': 0.013275560062690145, 'mean_word_length': 1.9991705833563727, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, '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': '4548000', 'n_tokens_mistral': 3568, 'n_tokens_neox': 3381, 'n_words': 905}