<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Visualize empirical correlation in HCP data
Step1:  Specify a folder containing preprocessed HCP data.
Step2:  Read data from file.
Step3:  This data has some negative values, which is unrealistic. This is taken care of by the function replaceNegativeData().
Step4:  Convert diffusion encoding to q-space features usable in prediction.
Step5:  Compute the spatial coordinates corresponding to each measurement.
Step6:  Compute a mask containing the brain. For this we first need to create an average b0-img, which is done by averaging all b0 images. Then, we pass this to createBrainMaskFromb0Data() which uses Dipy's median_otsu() function.
Step7:  Get the (x, y, z)-indices of voxels belonging to the mask using numpy.nonzero().
Step8:  Normalize by the b0-image, taking advantage of Numpy's broadcasting.
Step9:  Sort columns so that bvals are monotonically increasing
Step10:  Some measurements have b = 9950 or b = 10050. We want to treat them as part of the b=10000 shell.
Step11:  From the perspective of the qFeatures, each row in S is an observation and each column is a feature. Assuming that the rows are independent (enough) we can thus estimate the mean  signal as a function of qFeatures by averaging over the rows. The hypothesis is that the mean only depends on the q-magnitude, so we expect the empirical mean to be stepwise constant.
Step12:  The hypothesis seems to hold nicely. Now we can compute the empirical covariance.
Step13:  Investigate the covariance as a function of $\theta$ for each shell separately.
<ASSISTANT_TASK:>
Python Code:
import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import numpy as np
import matplotlib.pyplot as plt
from diGP.preprocessing import (readHCP,
                                averageb0Volumes,
                                createBrainMaskFromb0Data,
                                replaceNegativeData)
from diGP.dataManipulations import (log_q_squared,
                                    generateCoordinates)
%matplotlib inline
Explanation: Visualize empirical correlation in HCP data
End of explanation
dataPath = 'C:\\Users\\sesjojen\\Documents\\Data\\HumanConnectomeProject\\mgh_1007\\diff\\preproc'
print(dataPath)
Explanation: Specify a folder containing preprocessed HCP data.
End of explanation
gtab, data, voxelSize = readHCP(dataPath)
print(gtab.info)
Explanation: Read data from file.
End of explanation
data = replaceNegativeData(data, gtab)
Explanation: This data has some negative values, which is unrealistic. This is taken care of by the function replaceNegativeData().
End of explanation
def qMagnitudeTransform(q):
    return log_q_squared(q, c=1.)
qvecs = gtab.bvecs
qvals = gtab.qvals[:, np.newaxis]
qvalsTransformed = qMagnitudeTransform(qvals)
qFeatures = np.column_stack((qvalsTransformed, qvecs))
plt.hist(qvalsTransformed)
plt.xlabel('Transformed q-magnitude')
plt.show
Explanation: Convert diffusion encoding to q-space features usable in prediction.
End of explanation
spatialDimensions = data.shape[0:3]
coordinates = generateCoordinates(spatialDimensions, voxelSize)
coordinates.shape
Explanation: Compute the spatial coordinates corresponding to each measurement.
End of explanation
b0 = averageb0Volumes(data, gtab)
mask = createBrainMaskFromb0Data(b0)
zSlice = 60
f, (ax1, ax2) = plt.subplots(1, 2, sharey=True)
ax1.imshow(b0[:, :, zSlice])
ax2.imshow(mask[:, :, zSlice])
Explanation: Compute a mask containing the brain. For this we first need to create an average b0-img, which is done by averaging all b0 images. Then, we pass this to createBrainMaskFromb0Data() which uses Dipy's median_otsu() function.
End of explanation
maskIdx = np.nonzero(mask)
S = data[maskIdx[0], maskIdx[1], maskIdx[2], :]
Explanation: Get the (x, y, z)-indices of voxels belonging to the mask using numpy.nonzero().
End of explanation
S = S/b0[maskIdx[0], maskIdx[1], maskIdx[2], np.newaxis]
maxSignal = 1.5
percentExceedingMaxSignal = 100*np.sum(S>maxSignal)/np.prod(S.shape)
print('Replacing the top {} % values with {}.'.format(percentExceedingMaxSignal, maxSignal))
S[S > maxSignal] = maxSignal
Explanation: Normalize by the b0-image, taking advantage of Numpy's broadcasting.
End of explanation
bvalsIdxIncreasing = np.argsort(gtab.bvals)
bvals = gtab.bvals[bvalsIdxIncreasing]
bvecs = gtab.bvecs[bvalsIdxIncreasing, :]
S = S[:, bvalsIdxIncreasing]
qFeatures = qFeatures[bvalsIdxIncreasing, :]
Explanation: Sort columns so that bvals are monotonically increasing: 0, 1000, 3000, 5000, 10000. Adjust other data accordingly.
End of explanation
bvals[np.logical_and(bvals >= 9950., bvals <= 10050.)] = 10000. 
shells = np.unique(bvals)
print(shells)
shellIdx = [bvals == b for b in shells]
Explanation: Some measurements have b = 9950 or b = 10050. We want to treat them as part of the b=10000 shell.
End of explanation
mu = np.mean(S, axis=0)
plt.plot(mu)
Explanation: From the perspective of the qFeatures, each row in S is an observation and each column is a feature. Assuming that the rows are independent (enough) we can thus estimate the mean  signal as a function of qFeatures by averaging over the rows. The hypothesis is that the mean only depends on the q-magnitude, so we expect the empirical mean to be stepwise constant.
End of explanation
def computeEmpiricalCovariance(A):
    mu = np.mean(A, axis=0)
    Acentered = A - mu[np.newaxis, :]
    N = A.shape[0]
    empiricalCovariance = 1/(N-1)*np.dot(Acentered.T, Acentered)
    return empiricalCovariance
def removeShellMeans(A, shellIdx):
    A2 = np.zeros_like(A)
    for ind in shellIdx:
        shellMeans = np.mean(A[:, ind], axis=1)
        A2[:, ind] = A[:, ind] - shellMeans[:, np.newaxis]
    return A2
#R = np.corrcoef(removeShellMeans(S, shellIdx).T)
R = np.corrcoef(S.T)
plt.imshow(R, cmap='gray')
bvecsOuterProd = bvecs.dot(bvecs.T)
bvecsOuterProd[bvecsOuterProd > 1.] = 1.
bvecsOuterProd[bvecsOuterProd < -1.] = -1.
thetaMatrix = np.arccos(bvecsOuterProd)
thetaMatrix[thetaMatrix >= np.pi/2] -= np.pi; 
plt.imshow(thetaMatrix, cmap='gray')
Explanation: The hypothesis seems to hold nicely. Now we can compute the empirical covariance.
End of explanation
for i in np.arange(len(shells)):
    ind = np.squeeze(np.nonzero(shellIdx[i]))
    subMatrixIdx = np.ix_(ind, ind)
    ax = plt.subplot(2, 3, i+1)
    plt.scatter(thetaMatrix[subMatrixIdx]/np.pi, R[subMatrixIdx])
    plt.xlabel('$\\theta/\pi$')
    plt.ylabel('Correlation coefficient')
    plt.title('b = {}'.format(shells[i]))
    plt.xticks((-0.5, 0, 0.5))
    plt.tight_layout()
Explanation: Investigate the covariance as a function of $\theta$ for each shell separately.
End of explanation
<END_TASK>