<jupyter_start><jupyter_text>Master-M3R ComparisonIn this notebook, we will compare the individual step-by-step outputs of `master` and `m3r-release` to figure out where they diverge. For steps with a difference, we will additionally rerun each step with the previous matching inputs to figure out exactly where the two branches diverge.  PreprocessingIn this step, we will validate that the gradient tables are the same.<jupyter_code>m3r_grad_f = '/inputs/m3r-out/tmp/sub-0025427/ses-1/dwi/preproc/sub-0025427_ses-1_dwi_1.bvec'
m3r_t1_f = '/inputs/m3r-out/tmp/sub-0025427/ses-1/dwi/preproc/sub-0025427_ses-1_dwi_t1.nii.gz'

master_grad_f = '/inputs/ndmg-out/tmp/sub-0025427/ses-1/dwi/preproc/sub-0025427_ses-1_dwi_1.bvec'
master_t1_f = '/inputs/ndmg-out/tmp/sub-0025427/ses-1/dwi/preproc/sub-0025427_ses-1_dwi_t1.nii.gz'

# load the gradient tables
m3r_grad = np.loadtxt(m3r_grad_f, delimiter=' ')
master_grad = np.loadtxt(master_grad_f, delimiter=' ')
np.testing.assert_array_equal(m3r_grad, master_grad)

# load the t1ws
m3r_t1dat = nb.load(m3r_t1_f).get_data()
master_t1dat = nb.load(master_t1_f).get_data()
np.testing.assert_array_equal(m3r_t1dat, master_t1dat)<jupyter_output><empty_output><jupyter_text>RegistrationIn this step, we will validate that the registrations produce identical outputs, up to a rotation. Recall that the master branch uses the MNI brain not in MNI space; rather, the affine is flipped. To account for this, we will compare the outputs with and without an appropriate affine flip, to see whether FLIRT is robust to arbitrary rotations. If it is non-robust to rotations, we will use the master branch code coupled with the properly-oriented brains to compute the tractography and graphs to then further check these steps for equality between the two branches. Unrotated<jupyter_code>m3r_dwi_reg = '/inputs/m3r-out/dwi/registered/sub-0025427_ses-1_dwi_space-MNI152NLin6_res-1x1x1_registered.nii.gz'
master_dwi_reg = '/inputs/ndmg-out/reg/dwi/sub-0025427_ses-1_dwi_aligned.nii.gz'

m3r_dwidat = nb.load(m3r_dwi_reg).get_data()
master_dwidat = nb.load(master_dwi_reg).get_data()
np.testing.assert_array_equal(m3r_dwidat, master_dwidat)<jupyter_output><empty_output><jupyter_text>As we can see, this is because the affines are different:<jupyter_code>template = '/ndmg_atlases/atlas/MNI152NLin6_res-1x1x1_T1w.nii.gz'

print("MNI152 Affine")
print(nb.load(template).affine)
print("M3r Affine")
print(nb.load(m3r_dwi_reg).affine)
print("Master Affine")
print(nb.load(master_dwi_reg).affine)
ref = nb.load(template).get_data()

m3r2tempim = 'm3r2temp.png'
master2tempim = 'master2temp.png'

m3r2temp = plot_overlays(b0=ref, atlas=m3r_dwidat[:,:,:,0], edge=True)
plt.savefig(m3r2tempim)
plt.close()
master2temp = plot_overlays(b0=ref, atlas=master_dwidat[:,:,:,0], edge=True)
plt.savefig(master2tempim)
plt.close()

masterim = mpimg.imread(master2tempim)
m3rim = mpimg.imread(m3r2tempim)<jupyter_output><empty_output><jupyter_text>Below, we look at the diffusion brain's alignment with the MNI152 template at 1mm resolution and standard orientation:<jupyter_code>%matplotlib inline

dwi_cmp = side_by_side(masterim, m3rim, t1="Master to Temp", t2="M3r to Temp")
plt.close()
dwi_cmp<jupyter_output><empty_output><jupyter_text>As we can see, the master branch is backwards wrt the MNI152 template. Since the only difference between `master` and `m3r` through tractography steps is the registration being backwards on `master`, we continue using the tractography from `m3r`, noting that I did not touch anything for preprocessing, registration, and tractography on master with the exception of reorienting everything to standard space.<jupyter_code>tracks = np.load('/inputs/m3r-out/dwi/fiber/sub-0025427_ses-1_dwi_space-MNI152NLin6_res-1x1x1_fibers.npz')['arr_0']<jupyter_output><empty_output><jupyter_text>GraphFinally, we compare the graph generation code, using the graphs produced by `m3r` and comparing with the graphs produced from master's copy of the graph generation code. Here, I changed the following lines:```def make_graph(self, streamlines, attr=None):    """    Takes streamlines and produces a graph    **Positional Arguments:**            streamlines:                - Fiber streamlines either file or array in a dipy EuDX                  or compatible format.    """    nlines = np.shape(streamlines)[0]    print(" of Streamlines: " + str(nlines))    print_id = np.max((int(nlines*0.05), 1))   in case nlines*.05=0    for idx, streamline in enumerate(streamlines):        if (idx % print_id) == 0:            print(idx)        points = np.round(streamline).astype(int)        p = set()        for point in points:            try:                loc = self.rois[point[0], point[1], point[2]]            except IndexError:                pass            else:                pass            if loc:                p.add(loc)**      edges = set([tuple(sorted(x)) for x in product(p, p)])        for edge in edges:**          lst = tuple(sorted([str(node) for node in edge]))            self.edge_dict[lst] += 1**  edge_list = [(str(k[0]), str(k[1]), v) for k, v in self.edge_dict.items()]    self.g.add_weighted_edges_from(edge_list)```To:```def make_graph(self, streamlines, attr=None):    """    Takes streamlines and produces a graph    **Positional Arguments:**            streamlines:                - Fiber streamlines either file or array in a dipy EuDX                  or compatible format.    """    nlines = np.shape(streamlines)[0]    print(" of Streamlines: " + str(nlines))    for idx, streamline in enumerate(streamlines):        if (idx % int(nlines*0.05)) == 0:            print(idx)        points = np.round(streamline).astype(int)        p = set()        for point in points:            try:                loc = self.rois[point[0], point[1], point[2]]            except IndexError:                pass            else:                pass            if loc:                p.add(loc)**      edges = combinations(p, 2)        for edge in edges:**          lst = tuple([int(node) for node in edge])            self.edge_dict[tuple(sorted(lst))] += 1**  edge_list = [(k[0], k[1], v) for k, v in self.edge_dict.items()]    self.g.add_weighted_edges_from(edge_list)```Note the changed lines of interest are `**`'d.The below function is the original graphgen code from master:<jupyter_code>from __future__ import print_function

from itertools import combinations
from collections import defaultdict
import numpy as np
import networkx as nx
import nibabel as nb
import ndmg
import time


class graph_master(object):
    def __init__(self, N, rois, attr=None, sens="dwi"):
        """
        Initializes the graph with nodes corresponding to the number of ROIs
        **Positional Arguments:**
                N:
                    - Number of rois
                rois:
                    - Set of ROIs as either an array or niftii file)
                attr:
                    - Node or graph attributes. Can be a list. If 1 dimensional
                      will be interpretted as a graph attribute. If N
                      dimensional will be interpretted as node attributes. If
                      it is any other dimensional, it will be ignored.
        """
        self.N = N
        self.edge_dict = defaultdict(int)

        self.rois = nb.load(rois).get_data()
        n_ids = np.unique(self.rois)
        n_ids = n_ids[n_ids != 0]

        self.g = nx.Graph(name="Generated by NeuroData's MRI Graphs (ndmg)",
                          date=time.asctime(time.localtime()),
                          source="http://m2g.io",
                          region="brain",
                          sensor=sens,
                          ecount=0,
                          vcount=len(n_ids)
                          )
        print(self.g.graph)

        [self.g.add_node(ids) for ids in n_ids]
        pass

    def make_graph(self, streamlines, attr=None):
        """
        Takes streamlines and produces a graph
        **Positional Arguments:**
                streamlines:
                    - Fiber streamlines either file or array in a dipy EuDX
                      or compatible format.
        """
        nlines = np.shape(streamlines)[0]
        print("# of Streamlines: " + str(nlines))

        for idx, streamline in enumerate(streamlines):
            if (idx % int(nlines*0.05)) == 0:
                print(idx)

            points = np.round(streamline).astype(int)
            p = set()
            for point in points:
                try:
                    loc = self.rois[point[0], point[1], point[2]]
                except IndexError:
                    pass
                else:
                    pass

                if loc:
                    p.add(loc)

            edges = combinations(p, 2)
            for edge in edges:
                lst = tuple([int(node) for node in edge])
                self.edge_dict[tuple(sorted(lst))] += 1

        edge_list = [(k[0], k[1], v) for k, v in self.edge_dict.items()]
        self.g.add_weighted_edges_from(edge_list)

    def cor_graph(self, timeseries, attr=None):
        """
        Takes timeseries and produces a correlation matrix
        **Positional Arguments:**
            timeseries:
                -the timeseries file to extract correlation for.
                          dimensions are [numrois]x[numtimesteps]
        """
        print("Estimating correlation matrix for {} ROIs...".format(self.N))
        cor = np.corrcoef(timeseries)  # calculate pearson correlation

        roilist = np.unique(self.rois)
        roilist = roilist[roilist != 0]
        roilist = np.sort(roilist)

        for (idx_out, roi_out) in enumerate(roilist):
            for (idx_in, roi_in) in enumerate(roilist):
                self.edge_dict[tuple((roi_out, roi_in))] = float(np.absolute(
                    cor[idx_out, idx_in]))

        edge_list = [(k[0], k[1], v) for k, v in self.edge_dict.items()]

        self.g.add_weighted_edges_from(edge_list)
        pass

    def get_graph(self):
        """
        Returns the graph object created
        """
        try:
            return self.g
        except AttributeError:
            print("Error: the graph has not yet been defined.")
            pass

    def save_graph(self, graphname, fmt='edgelist'):
        """
        Saves the graph to disk
        **Positional Arguments:**
                graphname:
                    - Filename for the graph
        **Optional Arguments:**
                fmt:
                    - Output graph format
        """
        self.g.graph['ecount'] = nx.number_of_edges(self.g)
        g = nx.convert_node_labels_to_integers(self.g, first_label=1)
        if fmt == 'edgelist':
            nx.write_weighted_edgelist(g, graphname, encoding='utf-8')
        elif fmt == 'gpickle':
            nx.write_gpickle(g, graphname)
        elif fmt == 'graphml':
            nx.write_graphml(g, graphname)
        else:
            raise ValueError('edgelist, gpickle, and graphml currently supported')
        pass

    def summary(self):
        """
        User friendly wrapping and display of graph properties
        """
        print("\n Graph Summary:")
        print(nx.info(self.g))
        pass

parc = '/ndmg_atlases/label/desikan_res-1x1x1.nii.gz'
# make the graph with m3r's graphgen code, which is on branch m3r
m3r_graph = nd.graph(np.unique(nb.load(parc).get_data()) - 1, parc)
m3r_graph.make_graph(tracks)
# make the graph using master's original graphgen code above
master_graph = graph_master(np.unique(nb.load(parc).get_data()) - 1, parc)
master_graph.make_graph(tracks)
master_graph.g.nodes() == m3r_graph.g.nodes()
master_adj = nx.to_numpy_matrix(master_graph.g)
m3r_adj = nx.to_numpy_matrix(m3r_graph.g)
np.fill_diagonal(m3r_adj, 0)
np.testing.assert_equal(master_adj, m3r_adj)<jupyter_output><empty_output><jupyter_text>Showing that if the brains are registered properly on master (and hence, the tracts are properly registered) the resulting graphs are the same. Finally, just to be doubly sure, we will visualize the parcellation with MNI152:<jupyter_code>parc_dat = nb.load(parc).get_data()
f = plot_overlays(b0=ref, atlas=parc_dat, edge=True)
plt.close()
f<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10388127853881278, 'fraction_numerical': 0.020303326810176124, 'mean_word_length': 2.645957193816885, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, '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}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '584892', 'license': 'Apache-2.0', 'n_tokens_mistral': 4055, 'n_tokens_neox': 3850, 'n_words': 1015, 'path': 'notebooks/ebridge2/112/cmp_master_m3r.ipynb', 'repo_name': 'hhelm10/graspy', 'size': '12264'}