<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 In this notebook we will first implement a simplified version of the DCM model, in order to generate data for subsequent examples.
Step1:  Show the true graph
Step2:  Run an analysis using SPM
Now we will run an actual DCM analysis on these data using SPM.
1.Start spm using the command
Step3:  Click "done", which will take you to the next window to specify the B and C matrices.  Your entries should look like this
Step4:  This will create a file called "DCM_truemodel.mat" in the same directory as the other files.
Now estimate the model by clicking the DCM button, choosing "estimate" from the Action menu, and selecting "DCM_truemodel.mat" from the selector.  This will take a few minutes to compute.
Model 2
<ASSISTANT_TASK:>
Python Code:
import numpy
# use a consistent seed so that everyone has the same data
numpy.random.seed(1000)
import os,sys
%matplotlib inline
import matplotlib.pyplot as plt
sys.path.insert(0,'../utils')
from mkdesign import create_design_singlecondition
from graph_utils import show_graph_from_adjmtx
import math
from nipy.modalities.fmri.hemodynamic_models import spm_hrf,compute_regressor
import scipy.interpolate
results_dir = os.path.abspath("../results")
if not os.path.exists(results_dir):
    os.mkdir(results_dir)
# first let's build the model without the bilinear influence (aka PPI)
# after http://spm.martinpyka.de/?p=81
nregions=5
z=numpy.zeros(nregions)
# intrinsic connectivity
A=numpy.zeros((z.shape[0],z.shape[0]))
A=numpy.diag(numpy.ones(z.shape[0])*-1)
# add some structure
#A=A + numpy.diag(numpy.ones(z.shape[0]-1),k=-1)
A[2,1]=1
A[3,1]=1
B=numpy.zeros(A.shape)
B[2,0]=1
B[4,0]=1
C=numpy.zeros((z.shape[0],1))
C[0]=1
u=0
print (A)
print (B)
print (C)
# we are assuming a 1 second TR for the resulting data
# but the neural data are at a 1/16 millisecond time resolution
stepsize=.01
tslength=300
timepoints=numpy.arange(0,tslength,stepsize)
# create a blocked design
d,design=create_design_singlecondition(blockiness=1.0,deslength=tslength,blocklength=20,offset=20)
u=scipy.interpolate.griddata(numpy.arange(1,d.shape[0]),d,timepoints,fill_value=0)
def dcm_model(t,z,A,B,C,u):
    ut=numpy.abs(timepoints - t).argmin() 
    return (A.dot(z)+ u[ut]*B.dot(z) + C.dot(u[ut]).T)[0] 
def mk_dcm_dataset(timepoints,z,noise_sd):
    data=numpy.zeros((len(timepoints),len(z)))
    for i in range(1,len(timepoints)):
        data[i,:]=data[i-1,:] + dcm_model(timepoints[i],data[i-1,:],A,B,C,u)  + numpy.random.randn(len(z))*noise_sd 
    hrf=spm_hrf(stepsize,oversampling=1)
    data_conv=numpy.zeros(data.shape)
    for i in range(len(z)):
        data_conv[:,i]=numpy.convolve(data[:,i],hrf)[:data.shape[0]]        
    return data,data_conv    
noise_sd=2
data,data_conv=mk_dcm_dataset(timepoints,z,noise_sd)
numpy.savez(os.path.join(results_dir,'dcmdata.npz'),data=data_conv,A=A,B=B,C=C,u=u,d=d,design=design)
plt.subplot(211)
plt.plot(data_conv)
cc=numpy.corrcoef(data_conv.T)
print ('correlation matrix')
print (cc)
from sklearn.covariance import GraphLassoCV
import matplotlib.colors
glasso=GraphLassoCV()
glasso.fit(data_conv)
from pcor_from_precision import pcor_from_precision
pcor=pcor_from_precision(glasso.precision_)
print ('partial r^2 matrix')
print (pcor**2)
plt.figure(figsize=(10,5))
plt.subplot(141)
plt.imshow(A,interpolation='nearest',norm=matplotlib.colors.Normalize(vmin=-1,vmax=1))
plt.title('A mtx')
plt.subplot(142)
plt.imshow(B,interpolation='nearest',norm=matplotlib.colors.Normalize(vmin=-1,vmax=1))
plt.title('B mtx')
plt.subplot(143)
plt.imshow(cc,interpolation='nearest',norm=matplotlib.colors.Normalize(vmin=-1,vmax=1))
plt.title('correlation')
plt.subplot(144)
plt.imshow(pcor**2,interpolation='nearest',norm=matplotlib.colors.Normalize(vmin=-1,vmax=1))
plt.title('partial correlation')
Explanation: In this notebook we will first implement a simplified version of the DCM model, in order to generate data for subsequent examples.
End of explanation
gr=show_graph_from_adjmtx(A,B,C)
Explanation: Show the true graph
End of explanation
from IPython.display import Image
Image(filename='DCM_Amtx.png',retina=True) 
Explanation: Run an analysis using SPM
Now we will run an actual DCM analysis on these data using SPM.
1.Start spm using the command: 
spm fmri
2.Click the "Dynamic causal modelling" button on the main panel
3.From the "Action" menu, choose "specify".  
Model 1: True model
The selection menu should automatically go to the directory where the relevant files are (/home/vagrant/fmri-analysis-vm/analysis/connectivity/dcmfiles).  Assuming that it does, then first select the SPM.mat file listed in the directory, and click "Done" at the bottom.  The SPM.mat file contains a description of the experimental design that DCM needs.
You will be asked for the name: enter "truemodel"
Another file selector will appear, which is asking you to select the ROI files.  You should select all 5 of the VOI* files that are listed in the selector.  They need to be selected in order; the best way to do this is to hold the Shift key and click the last file in the list, which will select the entire set in order. Click "Done" again.  
You will be asked for a number of details about the design, for which you can just accept the default values.
Now we have to specify the three matrices in the SPM model (A,B, and C).  You will see a window that asks you to click radio buttons for the active connections.  For the A matrix, we need to include both the static and modulated connections from the A and B matrices above.  Your entry should look like this:
End of explanation
Image(filename='DCM_BCmtx.png',retina=True) 
Explanation: Click "done", which will take you to the next window to specify the B and C matrices.  Your entries should look like this:
End of explanation
Image(filename='DCM_BMS.png',retina=True) 
Explanation: This will create a file called "DCM_truemodel.mat" in the same directory as the other files.
Now estimate the model by clicking the DCM button, choosing "estimate" from the Action menu, and selecting "DCM_truemodel.mat" from the selector.  This will take a few minutes to compute.
Model 2: no modulation
Now let's specify the same model, but without any modulated connections.  We can do this using MATLAB code rather than doing it by hand. 
From the Utils menu, choose "Run M-file".  Using the file selector, go one directory up ('..') and select "mk_model_noppi.m". This will specify the model without modulation (i.e. with no nonzero entries in the B matrix) and then estimate it, saving the results to a file called "DCM_noppi.mat'.
Model 3: incorrect modulation
Now let's specify the same model, but with two incorrect modulated connections (so that it has the same number of parameters as the true model). Using the "Run M-file" command, select "mk_model_wrongppi.m"
Bayesian model comparison
A fundamental part of the DCM workflow is the comparison of models.  In this case, we can test for the PPI by comparing the three models that we just created.  
Choose "Dynamic causal modelling -> Action -> compare"
This will bring up a batch editor, in which we need to specify several variables.
Double click on "Directory" and select the current directory (".")
Click on "Inference method" and choose "Fixed effects" and click "specify"
Click "Data" then in the window below click "New subject". Then click ". Subject" under Data above, and choose "New Session" below.  Then click on "Models" above, click "Specify", and choose the the three DCM*.mat files that we created above, in the order truemodel, noppi, wrongppi.
Click the green arrow at the top to run the batch model.
End of explanation
<END_TASK>