Finished first pass of coherence analysis. Outputs now include arrays…

…, csv files and matplotlib-generated figures.

nipype/interfaces/nitime/analysis.py

@@ -11,6 +11,7 @@
 """
 
 import numpy as np
+import tempfile
 from nipype.utils.misc import package_check
 package_check('nitime')
 package_check('matplotlib')
@@ -19,10 +20,12 @@
 from nipype.interfaces.base import (TraitedSpec, File, InputMultiPath,
  OutputMultiPath, Undefined, traits,
  BaseInterface, isdefined)
+
+from nipype.utils.filemanip import fname_presuffix
+
 import nitime.analysis as nta
 from nitime.timeseries import TimeSeries
-
-from matplotlib.mlab import csv2rec
+import nitime.viz as viz
 
 class CoherenceAnalyzerInputSpec(TraitedSpec):
 
@@ -56,21 +59,26 @@ class CoherenceAnalyzerInputSpec(TraitedSpec):
  'which the analysis will average.',
  '[low,high] (Default [0.02,0.15]'))
 
- output_csv_file = File(desc='File to write output')
- output_figure_file = File(desc='File to write output figures')
- figure_type = traits.Enum('matrix','network',
+ output_csv_file = File(desc='File to write outputs (coherence,time-delay) with file-names: file_name_ {coherence,timedelay}')
+
+ output_figure_file = File(desc='File to write output figures (coherence,time-delay) with file-names: file_name_{coherence,timedelay}. Possible formats: .png,.svg,.pdf,.jpg,...')
+
+ figure_type = traits.Enum('matrix','network',usedefault=True,
  desc=("The type of plot to generate, where ",
  "'matrix' denotes a matrix image and",
- "'network' denotes a graph representation"))
+ "'network' denotes a graph representation.",
+ " Default: 'matrix'"))
 
 class CoherenceAnalyzerOutputSpec(TraitedSpec):
- coherence_array = traits.Array(desc=('The pairwise coherence values between',
- 'the ROIs'))
+ coherence_array = traits.Array(desc=('The pairwise coherence values',
+ 'between the ROIs'))
+
  timedelay_array = traits.Array(desc=('The pairwise time delays between the',
  'ROIs (in seconds)'))
 
  coherence_csv = File(desc = ('A csv file containing the pairwise ',
  'coherence values'))
+
  timedelay_csv = File(desc = ('A csv file containing the pairwise ',
  'time delay values'))
 
@@ -84,20 +92,25 @@ class CoherenceAnalyzer(BaseInterface):
  output_spec = CoherenceAnalyzerOutputSpec
 
  def _read_csv(self):
- """ Read from csv in_file and return a rec array """
+ """
+ Read from csv in_file and return an array and ROI names
 
-
+ The input file should have a first row containing the names of the
+ ROIs (strings)
 
+ the rest of the data will be read in and transposed so that the rows
+ (TRs) will becomes the second (and last) dimension of the array
+
+ """
  #Check that input conforms to expectations:
  first_row = open(self.inputs.in_file).readline()
  if not first_row[1].isalpha():
  raise ValueError("First row of in_file should contain ROI names as strings of characters")
 
-
  roi_names = open(self.inputs.in_file).readline().replace('\"','').strip('\n').split(',')
- data = np.loadtxt(self.inputs.in_file,skiprows=1,delimiter=',')
+ #Transpose, so that the time is the last dimension:
+ data = np.loadtxt(self.inputs.in_file,skiprows=1,delimiter=',').T
 
- #rec_array=csv2rec(self.inputs.in_file)
  return data,roi_names
 
  def _csv2ts(self):
@@ -124,10 +137,13 @@ def _run_interface(self, runtime):
  else:
  # get TS from inputs.in_TS
  TS = self.inputs.in_TS
- # deal with creating or storing ROI names
- if not TS.metadata.haskey('ROIs'):
- TS.metadata['ROIs']=['roi_%d' % x for x,_ in enumerate(TS.data)]
 
+ # deal with creating or storing ROI names:
+ if not TS.metadata.has_key('ROIs'):
+ self.ROIs=['roi_%d' % x for x,_ in enumerate(TS.data)]
+ else:
+ self.ROIs=TS.metadata['ROIs']
+
  A = nta.CoherenceAnalyzer(TS,
  method=dict(this_method='welch',
  NFFT=self.inputs.NFFT,
@@ -158,19 +174,74 @@ def _list_outputs(self):
  outputs['timedelay_array']=self.delay
 
  #Conditional
- if isdefined(self.inputs.output_csv_file):
- # we need to make a function that we call here that writes the coherence
- # values to this file "coherence_csv" and makes the time_delay csv file??
- outputs['coherence_csv'] = self.inputs.coherence_file
- outputs['timedelay_csv'] = timedelay_file
- if isdefined(self.inputs.output_figure_file):
- # we need to make a function that we call here that generates the
- # figure files for coherence and time delay (subplots, one fig??)
- # check self.inputs.figure_type but has a default??
- outputs['coherence_fig'] = self.inputs.output_figure_file
-
+ if isdefined(self.inputs.output_csv_file) and hasattr(self,'coherence'):
+ # we need to make a function that we call here that writes the
+ # coherence values to this file "coherence_csv" and makes the
+ # time_delay csv file??
+ self._make_output_files()
+ outputs['coherence_csv']=fname_presuffix(self.inputs.output_csv_file,suffix='_coherence')
+
+ outputs['timedelay_csv']=fname_presuffix(self.inputs.output_csv_file,suffix='_delay')
+
+ if isdefined(self.inputs.output_figure_file) and hasattr(self,
+ 'coherence'):
+ self._make_output_figures()
+ outputs['coherence_fig'] = fname_presuffix(self.inputs.output_figure_file,suffix='_coherence')
+ outputs['timedelay_fig'] = fname_presuffix(self.inputs.output_figure_file,suffix='_delay')
+
  return outputs
-
+ def _make_output_files(self):
+ """
+ Generate the output csv files.
+ """
+ for this in zip([self.coherence,self.delay],['coherence','delay']):
+ tmp_f = tempfile.mkstemp()[1]
+ np.savetxt(tmp_f,this[0],delimiter=',')
+
+ fid = open(fname_presuffix(self.inputs.output_csv_file,
+ suffix='_%s'%this[1]),'w+')
+ # this writes ROIs as header line
+ fid.write(','+','.join(self.ROIs)+'\n')
+ # this writes ROI and data to a line
+ for r, line in zip(self.ROIs, open(tmp_f)):
+ fid.write('%s,%s'%(r,line))
+ fid.close()
+
+
+ def _make_output_figures(self):
+ """
+ Generate the desired figure and save the files according to
+ self.inputs.output_figure_file
+
+ """
+
+ if self.inputs.figure_type == 'matrix':
+ fig_coh = viz.drawmatrix_channels(self.coherence,
+ channel_names=self.ROIs,
+ color_anchor=0)
+
+ fig_coh.savefig(fname_presuffix(self.inputs.output_figure_file,
+ suffix='_coherence'))
+
+ fig_dt = viz.drawmatrix_channels(self.delay,
+ channel_names=self.ROIs,
+ color_anchor=0)
+
+ fig_dt.savefig(fname_presuffix(self.inputs.output_figure_file,
+ suffix='_delay'))
+ else:
+ fig_coh = viz.drawgraph_channels(self.coherence,
+ channel_names=self.ROIs)
+
+ fig_coh.savefig(fname_presuffix(self.inputs.output_figure_file,
+ suffix='_coherence'))
+
+ fig_dt = viz.drawgraph_channels(self.delay,
+ channel_names=self.ROIs)
+
+ fig_dt.savefig(fname_presuffix(self.inputs.output_figure_file,
+ suffix='_delay'))
+
 class GetTimeSeriesInputSpec():
  pass
 class GetTimeSeriesOutputSpec():

nipype/interfaces/nitime/tests/test_nitime.py

@@ -4,13 +4,20 @@
 import tempfile
 import shutil
 
+import numpy as np
+
 from nipype.testing import (assert_equal, assert_not_equal, assert_raises,
  with_setup, TraitError, parametric, skipif)
 
 from nipype.testing import example_data
 
 import nipype.interfaces.nitime as nitime 
 
+import nitime.analysis as nta
+import nitime.timeseries as ts
+
+from matplotlib.mlab import csv2rec
+
 def test_read_csv():
  """Test that reading the data from csv file gives you back a reasonable
  time-series object """
@@ -24,3 +31,42 @@ def test_read_csv():
  data,roi_names = CA._read_csv()
  yield assert_equal, data[0][0],10125.9
  yield assert_equal, roi_names[0],'WM' 
+
+def test_coherence_analysis():
+ """Test that the coherence analyzer works """
+
+ #This is the nipype interface analysis:
+ CA = nitime.CoherenceAnalyzer()
+ CA.inputs.TR = 1.89
+ CA.inputs.in_file = example_data('fmri_timeseries.csv')
+ o = CA.run()
+ yield assert_equal,o.outputs.coherence_array.shape,(31,31)
+
+ #This is the nitime analysis:
+ TR=1.89
+ data_rec = csv2rec(example_data('fmri_timeseries.csv'))
+ roi_names= np.array(data_rec.dtype.names)
+ n_samples = data_rec.shape[0]
+ data = np.zeros((len(roi_names),n_samples))
+
+ for n_idx, roi in enumerate(roi_names):
+ data[n_idx] = data_rec[roi]
+
+ T = ts.TimeSeries(data,sampling_interval=TR)
+
+ yield assert_equal,CA._csv2ts().data,T.data
+
+ T.metadata['roi'] = roi_names 
+ C = nta.CoherenceAnalyzer(T,method=dict(this_method='welch',
+ NFFT=CA.inputs.NFFT,
+ n_overlap=CA.inputs.n_overlap))
+
+ freq_idx = np.where((C.frequencies>CA.inputs.frequency_range[0]) *
+ (C.frequencies<CA.inputs.frequency_range[1]))[0]
+
+ #Extract the coherence and average across these frequency bands: 
+ coh = np.mean(C.coherence[:,:,freq_idx],-1) #Averaging on the last dimension 
+
+ yield assert_equal,o.outputs.coherence_array,coh
+
+