GitHub PullRequestReviewCommentEvent at 2014-02-19 13:39:20
Repository: nilearn
Actor: VirgileFritsch
Payload: {"comment":{"url":"https://api.github.com/repos/nilearn/nilearn/pulls/comments/9864572","id":9864572,"diff_hunk":"@@ -0,0 +1,593 @@\n+\"\"\"\n+Massively Univariate Linear Model estimated with OLS and permutation test.\n+\n+\"\"\"\n+# Author: Benoit Da Mota, <benoit.da_mota@inria.fr>, sept. 2011\n+# refactorized by Virgile Fritsch, <virgile.fritsch@inria.fr>, jan. 2014\n+import warnings\n+import numpy as np\n+from scipy import linalg, stats\n+from sklearn.utils import gen_even_slices, check_random_state\n+import sklearn.externals.joblib as joblib\n+\n+\n+def normalize_matrix_on_axis(m, axis=0):\n+    \"\"\" Normalize a 2D matrix on an axis.\n+\n+    Parameters\n+    ----------\n+    m : numpy 2D array,\n+        The matrix to normalize\n+    axis : integer in {0, 1}, optional\n+        A valid axis to normalize accross (0 by default)\n+\n+    Returns\n+    -------\n+    ret : numpy array, shape = m.shape\n+        The normalize matrix\n+\n+    Examples\n+    --------\n+    >>> import numpy as np\n+    >>> from nilearn.mass_univariate.permuted_least_squares import (\n+    ...     normalize_matrix_on_axis)\n+    >>> X = np.array([[0, 4], [1, 0]])\n+    >>> normalize_matrix_on_axis(X)\n+    array([[ 0.,  1.],\n+           [ 1.,  0.]])\n+    >>> normalize_matrix_on_axis(X, axis=1)\n+    array([[ 0.,  1.],\n+           [ 1.,  0.]])\n+\n+    \"\"\"\n+    if m.ndim > 2:\n+        raise ValueError('This function only accepts 2D arrays. '\n+                         'An array of shape %r was passed.' % m.shape)\n+\n+    if axis == 0:\n+        ret = m / np.sqrt(np.sum(m ** 2, axis=0))\n+    elif axis == 1:\n+        ret = normalize_matrix_on_axis(m.T).T\n+    else:\n+        raise ValueError('Invalid axis in normalization: %r' % axis)\n+    return np.ascontiguousarray(ret)\n+\n+\n+def orthonormalize_matrix(m, tol=1.e-12):\n+    \"\"\" Orthonormalize a matrix and complete it with zeros to preserve shape.\n+\n+    Parameters\n+    ----------\n+    m : numpy array,\n+      The matrix to orthonormalize\n+\n+    Returns\n+    -------\n+    ret : numpy array, shape = m.shape\n+      The orthonormalized matrix. It may be completed with zeros in order\n+      to preserve the initial matrix shape.\n+\n+    Examples\n+    --------\n+    >>> import numpy as np\n+    >>> from nilearn.mass_univariate.permuted_least_squares import (\n+    ...     orthonormalize_matrix)\n+    >>> X = np.array([[0, 1], [0, 1]])\n+    >>> orthonormalize_matrix(X)\n+    array([[ 0.70710678,  0.        ],\n+           [ 0.70710678,  0.        ]])\n+    >>> X = np.array([[0, 1], [4, 0]])\n+    >>> orthonormalize_matrix(X)\n+    array([[ 0., -1.],\n+           [-1.,  0.]])\n+\n+    \"\"\"\n+    U, s, _ = linalg.svd(m, full_matrices=False)\n+    n_eig = s[abs(s) > tol].size\n+    tmp = np.dot(U, np.diag(s))[:, :n_eig]\n+    n_null_eig = s.size - n_eig\n+    tmp = normalize_matrix_on_axis(tmp)\n+    if n_null_eig > 0:\n+        tmp = np.ascontiguousarray(\n+            np.hstack((tmp, np.zeros((tmp.shape[0], n_null_eig)))))\n+    return tmp\n+\n+\n+class GrowableSparseArray(object):\n+    \"\"\"Data structure to contain data from numerous estimations.\n+\n+    Examples of application are all resampling schemes\n+    (bootstrap, permutations, ...)\n+\n+    GrowableSparseArray can be seen as a three-dimensional array that contains\n+    scores associated with three position indices corresponding to\n+    (i) an iteration (or estimation), (ii) a test variate and\n+    (iii) a target variate.\n+    Memory is pre-allocated to store a large number of scores. The structure\n+    can be indexed efficiently according to three dimensions to add new\n+    scores at the right position fast.\n+    The allocated space can be extended if needed, but we want to avoid this\n+    because it is costly. User should carefully initialize the structure.\n+    Only scores above a predetermined threshold are actually stored, others\n+    are ignored.\n+\n+    Attributes\n+    ----------\n+    n_elts: int\n+      The total number of scores actually stored into the data structure\n+    n_iter: int\n+      Number of trials (using as many iterators)\n+    max_elts: int\n+      Maximum number of scores that can be stored into the structure\n+    data: array-like, own-designed dtype\n+      The actual scores corresponding to all the estimations.\n+      dtype is built so that every score is associated with three position\n+      GrowableSparseArray can be seen as a three-dimensional array so every\n+      score is associated with three position indices corresponding to\n+      (i) an iteration (or an estimator) ('iter_id'),\n+      (ii) a test variate ('x_id') and\n+      (iii) a target variate ('y_id').\n+    sizes: array-like, shape=(n_iter, )\n+      The number of scores stored for each estimation.\n+      Useful to select a range of values from iteration ids.\n+    threshold: float,\n+      Sparsity threshold used to discard scores that are to low to have a\n+      chance to correspond to a maximum value amongst all the scores of\n+      a given iteration.\n+\n+    \"\"\"\n+    def __init__(self, n_iter=10000, n_elts=0, max_elts=None,\n+                 threshold=-np.inf):\n+        self.n_elts = n_elts\n+        self.n_iter = n_iter\n+        self.max_elts = max(max_elts, n_elts)\n+        self.data = np.empty(\n+            self.max_elts,\n+            dtype=[('iter_id', np.int32), ('x_id', np.int32),\n+                   ('y_id', np.int32), ('score', np.float32)])\n+        self.sizes = np.zeros(n_iter, dtype=int)\n+        self.threshold = threshold\n+\n+    def get_data(self):\n+        return self.data[:self.n_elts]\n+\n+    def merge(self, others):\n+        \"\"\"Copy one or several GrowableSparseArray into the current structure.\n+\n+        Parameters\n+        ----------\n+        others: list of GrowableSparseArray or GrowableSparseArray\n+          The structures to be merged into the current structure.\n+\n+        \"\"\"\n+        if isinstance(others, GrowableSparseArray):\n+            return self.merge([others])\n+        if not isinstance(others, list) and not isinstance(others, tuple):\n+            raise TypeError(\n+                '\\'others\\' is not a list/tuple of GrowableSparseArray '\n+                'or a GrowableSparseArray.')\n+        for gs_array in others:\n+            if not isinstance(gs_array, GrowableSparseArray):\n+                raise TypeError('List element is not a GrowableSparseArray.')\n+            if gs_array.n_iter != self.n_iter:\n+                raise ValueError('Cannot merge a structure with %d iterations '\n+                                'into a structure with %d iterations.'\n+                                % (gs_array.n_iter, self.n_iter))\n+\n+        acc_sizes = [self.sizes]\n+        acc_data = [self.get_data()]\n+        for gs_array in others:\n+            # threshold the data to respect self.threshold\n+            if gs_array.threshold < self.threshold:\n+                gs_array_data_thresholded = (\n+                    gs_array.get_data()[gs_array.get_data()['score']\n+                                       >= self.threshold])\n+                acc_sizes.append([gs_array_data_thresholded.size])\n+                acc_data.append(gs_array_data_thresholded)\n+            elif gs_array.threshold > self.threshold:\n+                warnings.warn('Merging a GrowableSparseArray into another '\n+                              'with a lower threshold: parent array may '\n+                              'contain less scores than its threshold '\n+                              'suggests.')\n+                acc_sizes.append(gs_array.sizes)\n+                acc_data.append(gs_array.get_data())\n+            else:\n+                acc_sizes.append(gs_array.sizes)\n+                acc_data.append(gs_array.get_data())\n+\n+        self.sizes = np.array(acc_sizes).sum(axis=0)\n+        self.data = np.concatenate(acc_data)\n+        self.n_elts = self.sizes.sum()\n+        self.max_elts = self.n_elts\n+        self.data = np.sort(self.data, order=['iter_id', 'x_id', 'y_id'])\n+\n+    def append(self, iter_id, iter_data, y_offset=0):\n+        \"\"\"Add the data of one estimation (iteration) into the structure.\n+\n+        This is done in a memory-efficient way, by taking into account\n+        pre-allocated space.\n+\n+        Parameters\n+        ----------\n+        iter_id: int,\n+          ID of the estimation we are inserting into the structure\n+        iter_data: array-like, shape=(n_targets_chunk, n_regressors)\n+          Scores corresponding to the iteration chunk to be inserted into\n+          the data structure.\n+        y_offset: int,\n+          Position of the target variates chunk relative to the original\n+          dataset.\n+\n+        \"\"\"\n+        # we only store float32 to save space\n+        iter_data = iter_data.astype('float32')\n+        # we sparsify the matrix wrt. threshold using coordinates list\n+        y_idx, x_idx = (iter_data >= self.threshold).nonzero()\n+        score_size = len(x_idx)\n+        new_n_elts = score_size + self.n_elts\n+        if (new_n_elts > self.max_elts or\n+            self.sizes[iter_id + 1:].sum() > 0):  # insertion (costly)\n+            new_data = np.empty(score_size,\n+                        dtype=[('iter_id', np.int32), ('x_id', np.int32),\n+                               ('y_id', np.int32), ('score', np.float32)])\n+            new_data['x_id'][:] = x_idx\n+            new_data['y_id'][:] = y_idx + y_offset\n+            new_data['score'][:] = iter_data[y_idx, x_idx]\n+            new_data['iter_id'][:] = iter_id\n+            gs_array = GrowableSparseArray(\n+                self.n_iter, threshold=self.threshold)\n+            gs_array.data = new_data\n+            gs_array.sizes = np.zeros((gs_array.n_iter), dtype=int)\n+            gs_array.sizes[iter_id] = score_size\n+            gs_array.n_elts = score_size\n+            gs_array.max_elts = score_size\n+            self.merge(gs_array)\n+        else:  # it fits --> updates (efficient)\n+            self.data['x_id'][self.n_elts:new_n_elts] = x_idx\n+            self.data['y_id'][self.n_elts:new_n_elts] = y_idx + y_offset\n+            self.data['score'][self.n_elts:new_n_elts] = (\n+                iter_data[y_idx, x_idx])\n+            self.data['iter_id'][self.n_elts:new_n_elts] = iter_id\n+            self.sizes[iter_id] += score_size\n+            self.n_elts = new_n_elts\n+\n+\n+def f_score(vars1, vars2, covars, lost_dof):\n+    \"\"\"Compute F-score associated with the regression of vars2 against vars1\n+\n+    Covariates are taken into account\n+\n+    Parameters\n+    ----------\n+    vars1: array-like, shape=(n_samples, n_var1)\n+      Explanatory variates\n+    vars2: array-like, shape=(n_var2, n_samples)\n+      Targets variates\n+    covars, array-like, shape=(n_samples, n_covars) or None\n+      Confounding variates\n+    lost_dof: int,\n+      Lost degrees of freedom\n+\n+    Returns\n+    -------\n+    score: array-like, shape=(n_var2, n_var1)\n+      F-scores associated with the tests of each explanatory variate against\n+      each target variate (in the presence of covars).\n+\n+    \"\"\"\n+    dof = vars2.shape[1] - 1 - lost_dof\n+    beta_vars2_vars1 = np.dot(vars2, vars1)\n+    b2 = beta_vars2_vars1 ** 2\n+    if covars is None:\n+        rss = (1 - b2)\n+    else:\n+        beta_vars2_covars = np.dot(vars2, covars)\n+        a2 = np.sum(beta_vars2_covars ** 2, 1)\n+        rss = (1 - a2[:, np.newaxis] - b2)\n+    score = b2 / rss\n+    score *= dof\n+    return score\n+\n+\n+def _permuted_ols_on_chunk(tested_vars, target_vars_chunk,\n+                           confounding_vars=None, n_perm=10000,\n+                           sparsity_threshold=1e-04,\n+                           target_vars_chunk_position=0,\n+                           intercept_test=True, random_state=0):","path":"nilearn/mass_univariate/permuted_least_squares.py","position":296,"original_position":296,"commit_id":"3fefed06dff95a8d21508a634f72394f1e148c18","original_commit_id":"3fefed06dff95a8d21508a634f72394f1e148c18","user":{"login":"VirgileFritsch","id":263280,"avatar_url":"https://gravatar.com/avatar/06ae2d71d426285c4cf2206301581bb0?d=https%3A%2F%2Fidenticons.github.com%2F9e896193cb8d51362144b9a7f427f64b.png&r=x","gravatar_id":"06ae2d71d426285c4cf2206301581bb0","url":"https://api.github.com/users/VirgileFritsch","html_url":"https://github.com/VirgileFritsch","followers_url":"https://api.github.com/users/VirgileFritsch/followers","following_url":"https://api.github.com/users/VirgileFritsch/following{/other_user}","gists_url":"https://api.github.com/users/VirgileFritsch/gists{/gist_id}","starred_url":"https://api.github.com/users/VirgileFritsch/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/VirgileFritsch/subscriptions","organizations_url":"https://api.github.com/users/VirgileFritsch/orgs","repos_url":"https://api.github.com/users/VirgileFritsch/repos","events_url":"https://api.github.com/users/VirgileFritsch/events{/privacy}","received_events_url":"https://api.github.com/users/VirgileFritsch/received_events","type":"User","site_admin":false},"body":"Please stop making me look like a pretentious newcomer that want to impress by its skills and self-pretending revolutionary ideas. If this is what you think, please immediately close this PR so we will stop trying to make quality code together.\r\nYou seem so attached to recognizing your favorite conventions, and so eager to find it in my code that it completely blurred your understanding of what my code was actually doing.\r\nAll I wish is to respect conventions the best I can and I am not responsible if you spend your day fighting people that do not have this wish. So please remind that we are working in the same direction.\r\n\r\nI am trying to tell you that `random_seed` has been set up to control the state of a routine which is supposed to rely on random numbers, especially for the sake of debugging, testing, and examples generation.\r\nHere, the algorithm *requires* that we *always* have the control of the random numbers because they should be the same in simultaneous calls to the routine that uses them. It is not something I do for my pleasure, it is permutation tests theory that states that the F-max correction implies all the targets are permuted in the exact same fashion to enable comparison between them (finding the max).\r\n\r\nI can put `random_state=None` by default, but then I would immediately assign it a fixed value in the first function instruction (yielding a fake impression to the code reviewer).\r\n\r\nSo:\r\n> why? [renaming `random_state` to `seed`]\r\n\r\nTo avoid conventional people --that recognize a familiar name-- think that I did not follow the conventions whereas I was actually doing something different.\r\n\r\n> what is the gain?\r\n\r\nIt becomes a compulsory argument, so that nobody uses the code without controlling the seed, and avoid making the mistake which whould consist of permuting chunks differently.\r\n\r\nThe big picture of what the permuted_ols code would become:\r\n1. `permuted_ols` has a `random_state` argument that is defaulted to `None`\r\n2. ```rng = check_random_state(random_state)```\r\n3. ```seed = rng.randint(sys.maxint)  # the common seed itself is random, or controlled by `random_state` in examples and tests```\r\n4. parallel calls to `_permuted_ols_on_chunk(seed, ...)`  (all chunks have the same seed, they cannot have a different one)\r\n\r\nNow, if you have a better option...","created_at":"2014-02-19T13:39:20Z","updated_at":"2014-02-19T13:39:20Z","html_url":"https://github.com/nilearn/nilearn/pull/157#discussion_r9864572","pull_request_url":"https://api.github.com/repos/nilearn/nilearn/pulls/157","_links":{"self":{"href":"https://api.github.com/repos/nilearn/nilearn/pulls/comments/9864572"},"html":{"href":"https://github.com/nilearn/nilearn/pull/157#discussion_r9864572"},"pull_request":{"href":"https://api.github.com/repos/nilearn/nilearn/pulls/157"}}}}