topo.analysis.vision

1 """ 2 Vision-specific analysis functions. 3 4 $Id: featureresponses.py 7714 2008-01-24 16:42:21Z antolikjan $ 5 """ 6 __version__='$Revision: 7714 $' 7 8 from math import pi,sin,cos 9 10 import numpy 11 from numpy.oldnumeric import Float 12 from numpy import zeros, array, size, object_ 13 #import scipy 14 15 import param 16 from param import normalize_path 17 18 try: 19 import matplotlib 20 import pylab 21 except ImportError: 22 param.Parameterized(name=__name__).warning("Could not import matplotlib; module will not be useable.") 23 from topo.command.basic import ImportErrorRaisingFakeModule 24 pylab = ImportErrorRaisingFakeModule("matplotlib") 25 26 # CEBALERT: commands in here should inherit from the appropriate base 27 # class (Command or PylabPlotCommand). 28 29 import topo 30 from topo.base.cf import CFSheet 31 from topo.base.sheetview import SheetView 32 from topo.plotting.plotgroup import create_plotgroup 33 from topo.command.analysis import measure_sine_pref 34 35 from topo import numbergen 36 37 max_value = 0 38 global_index = () 39

40 -def _complexity_rec(x,y,index,depth,fm):

41 """ 42 Recurrent helper function for complexity() 43 """ 44 global max_value 45 global global_index 46 if depth<size(fm.features): 47 for i in range(size(fm.features[depth].values)): 48 _complexity_rec(x,y,index + (i,),depth+1,fm) 49 else: 50 if max_value < fm.full_matrix[index][x][y]: 51 global_index = index 52 max_value = fm.full_matrix[index][x][y]

53 54 55

56 -def complexity(full_matrix):

57 global global_index 58 global max_value 59 """This function expects as an input a object of type FullMatrix which contains 60 responses of all neurons in a sheet to stimuly with different varying parameter values. 61 One of these parameters (features) has to be phase. In such case it computes the classic 62 modulation ratio (see Hawken et al. for definition) for each neuron and returns them as a matrix. 63 """ 64 rows,cols = full_matrix.matrix_shape 65 complexity = zeros(full_matrix.matrix_shape) 66 complex_matrix = zeros(full_matrix.matrix_shape,object_) 67 fftmeasure = zeros(full_matrix.matrix_shape,Float) 68 i = 0 69 for f in full_matrix.features: 70 if f.name == "phase": 71 72 phase_index = i 73 break 74 i=i+1 75 sum = 0.0 76 res = 0.0 77 average = 0.0 78 for x in range(rows): 79 for y in range(cols): 80 complex_matrix[x,y] = []# 81 max_value=-0.01 82 global_index = () 83 _complexity_rec(x,y,(),0,full_matrix) 84 85 #compute the sum of the responses over phases given the found index of highest response 86 87 iindex = array(global_index) 88 sum = 0.0 89 for i in range(size(full_matrix.features[phase_index].values)): 90 iindex[phase_index] = i 91 sum = sum + full_matrix.full_matrix[tuple(iindex.tolist())][x][y] 92 93 #average 94 average = sum / float(size(full_matrix.features[phase_index].values)) 95 96 res = 0.0 97 #compute the sum of absolute values of the responses minus average 98 for i in range(size(full_matrix.features[phase_index].values)): 99 iindex[phase_index] = i 100 res = res + abs(full_matrix.full_matrix[tuple(iindex.tolist())][x][y] - average) 101 complex_matrix[x,y] = complex_matrix[x,y] + [full_matrix.full_matrix[tuple(iindex.tolist())][x][y]] 102 103 if x==43 and y==43: 104 pylab.figure() 105 ax = pylab.subplot(111) 106 z = complex_matrix[x,y][:] 107 z.append(z[0]) 108 pylab.plot(z,linewidth=4) 109 pylab.axis(xmin=0.0,xmax=numpy.pi) 110 ax.yaxis.set_major_locator(matplotlib.ticker.MaxNLocator(4)) 111 pylab.xticks([0,len(z)/2,len(z)-1], ['0','pi/2','pi']) 112 pylab.savefig(normalize_path(str(topo.sim.time()) + str(complex_matrix[x,y][0])+ 'modulation_response[43,43].png')) 113 114 if x==45 and y==45: 115 pylab.figure() 116 ax = pylab.subplot(111) 117 z = complex_matrix[x,y][:] 118 z.append(z[0]) 119 pylab.plot(z,linewidth=4) 120 pylab.axis(xmin=0.0,xmax=numpy.pi) 121 ax.yaxis.set_major_locator(matplotlib.ticker.MaxNLocator(4)) 122 pylab.xticks([0,len(z)/2,len(z)-1], ['0','pi/2','pi']) 123 pylab.savefig(normalize_path(str(topo.sim.time()) + str(complex_matrix[x,y][0])+ 'modulation_response[45,45].png')) 124 125 fft = numpy.fft.fft(complex_matrix[x,y]+complex_matrix[x,y]+complex_matrix[x,y]+complex_matrix[x,y],2048) 126 first_har = 2048/len(complex_matrix[0,0]) 127 if abs(fft[0]) != 0: 128 fftmeasure[x,y] = 2 *abs(fft[first_har]) /abs(fft[0]) 129 else: 130 fftmeasure[x,y] = 0 131 return fftmeasure

132 133

134 -def compute_ACDC_orientation_tuning_curves(full_matrix,curve_label,sheet):

135 136 """ This function allows and alternative computation of orientation tuning curve where 137 for each given orientation the response is computed as a maximum of AC or DC component 138 across the phases instead of the maximum used as a standard in Topographica""" 139 # this method assumes that only single frequency has been used 140 i = 0 141 for f in full_matrix.features: 142 if f.name == "phase": 143 phase_index = i 144 if f.name == "orientation": 145 orientation_index = i 146 if f.name == "frequency": 147 frequency_index = i 148 i=i+1 149 print sheet.curve_dict 150 if not sheet.curve_dict.has_key("orientationACDC"): 151 sheet.curve_dict["orientationACDC"]={} 152 sheet.curve_dict["orientationACDC"][curve_label]={} 153 154 rows,cols = full_matrix.matrix_shape 155 for o in xrange(size(full_matrix.features[orientation_index].values)): 156 s_w = zeros(full_matrix.matrix_shape) 157 for x in range(rows): 158 for y in range(cols): 159 or_response=[] 160 for p in xrange(size(full_matrix.features[phase_index].values)): 161 index = [0,0,0] 162 index[phase_index] = p 163 index[orientation_index] = o 164 index[frequency_index] = 0 165 or_response.append(full_matrix.full_matrix[tuple(index)][x][y]) 166 167 fft = numpy.fft.fft(or_response+or_response+or_response+or_response,2048) 168 first_har = 2048/len(or_response) 169 s_w[x][y] = numpy.maximum(2 *abs(fft[first_har]),abs(fft[0])) 170 s = SheetView((s_w,sheet.bounds), sheet.name , sheet.precedence, topo.sim.time(),sheet.row_precedence) 171 sheet.curve_dict["orientationACDC"][curve_label].update({full_matrix.features[orientation_index].values[o]:s})

172 173 174

175 -def phase_preference_scatter_plot(sheet_name,diameter=0.39):

176 r = numbergen.UniformRandom(seed=1023) 177 preference_map = topo.sim[sheet_name].sheet_views['PhasePreference'] 178 offset_magnitude = 0.03 179 datax = [] 180 datay = [] 181 (v,bb) = preference_map.view() 182 for z in zeros(66): 183 x = (r() - 0.5)*2*diameter 184 y = (r() - 0.5)*2*diameter 185 rand = r() 186 xoff = sin(rand*2*pi)*offset_magnitude 187 yoff = cos(rand*2*pi)*offset_magnitude 188 xx = max(min(x+xoff,diameter),-diameter) 189 yy = max(min(y+yoff,diameter),-diameter) 190 x = max(min(x,diameter),-diameter) 191 y = max(min(y,diameter),-diameter) 192 [xc1,yc1] = topo.sim[sheet_name].sheet2matrixidx(xx,yy) 193 [xc2,yc2] = topo.sim[sheet_name].sheet2matrixidx(x,y) 194 if((xc1==xc2) & (yc1==yc2)): continue 195 datax = datax + [v[xc1,yc1]] 196 datay = datay + [v[xc2,yc2]] 197 198 for i in range(0,len(datax)): 199 datax[i] = datax[i] * 360 200 datay[i] = datay[i] * 360 201 if(datay[i] > datax[i] + 180): datay[i]= datay[i]- 360 202 if((datax[i] > 180) & (datay[i]> 180)): datax[i] = datax[i] - 360; datay[i] = datay[i] - 360 203 if((datax[i] > 180) & (datay[i] < (datax[i]-180))): datax[i] = datax[i] - 360; #datay[i] = datay[i] - 360 204 205 f = pylab.figure() 206 ax = f.add_subplot(111, aspect='equal') 207 pylab.plot(datax,datay,'ro') 208 pylab.plot([0,360],[-180,180]) 209 pylab.plot([-180,180],[0,360]) 210 pylab.plot([-180,-180],[360,360]) 211 ax.axis([-180,360,-180,360]) 212 pylab.xticks([-180,0,180,360], [-180,0,180,360]) 213 pylab.yticks([-180,0,180,360], [-180,0,180,360]) 214 pylab.grid() 215 pylab.savefig(normalize_path(str(topo.sim.timestr()) + sheet_name + "_scatter.png"))

216 217 218 219 ############################################################################### 220 # JABALERT: Should we move this plot and command to analysis.py or 221 # pylabplot.py, where all the rest are? 222 # 223 # In any case, it requires generalization; it should not be hardcoded 224 # to any particular map name, and should just do the right thing for 225 # most networks for which it makes sense. E.g. it already measures 226 # the ComplexSelectivity for all measured_sheets, but then 227 # plot_modulation_ratio only accepts two with specific names. 228 # plot_modulation_ratio should just plot whatever it is given, and 229 # then analyze_complexity can simply pass in whatever was measured, 230 # with the user controlling what is measured using the measure_map 231 # attribute of each Sheet. That way the complexity of any sheet could 232 # be measured, which is what we want. 233 # 234 # Specific changes needed: 235 # - Make plot_modulation_ratio accept a list of sheets and 236 # plot their individual modulation ratios and combined ratio. 237 # - Remove complex_sheet_name argument, which is no longer needed 238 # - Make sure it still works fine even if V1Simple doesn't exist; 239 # as this is just for an optional scatter plot, it's fine to skip 240 # it. 241 # - Preferably remove the filename argument by default, so that 242 # plots will show up in the GUI 243 244

245 -def analyze_complexity(full_matrix,simple_sheet_name,complex_sheet_name,filename=None):

246 """ 247 Compute modulation ratio for each neuron, to distinguish complex from simple cells. 248 249 Uses full_matrix data obtained from measure_or_pref(). 250 251 If there is a sheet named as specified in simple_sheet_name, 252 also plots its phase preference as a scatter plot. 253 """ 254 import topo 255 measured_sheets = [s for s in topo.sim.objects(CFSheet).values() 256 if hasattr(s,'measure_maps') and s.measure_maps] 257 258 for sheet in measured_sheets: 259 # Divide by two to get into 0-1 scale - that means simple/complex boundry is now at 0.5 260 complx = array(complexity(full_matrix[sheet]))/2.0 261 # Should this be renamed to ModulationRatio? 262 sheet.sheet_views['ComplexSelectivity']=SheetView((complx,sheet.bounds), sheet.name , sheet.precedence, topo.sim.time(),sheet.row_precedence) 263 import topo.command.pylabplot 264 topo.command.pylabplot.plot_modulation_ratio(full_matrix,simple_sheet_name=simple_sheet_name,complex_sheet_name=complex_sheet_name,filename=filename) 265 266 # Avoid error if no simple sheet exists 267 try: 268 phase_preference_scatter_plot(simple_sheet_name,diameter=0.24999) 269 except AttributeError: 270 print "Skipping phase preference scatter plot; could not analyze region %s." \ 271 % simple_sheet_name

272 273

274 -class measure_and_analyze_complexity(measure_sine_pref):

275 """Macro for measuring orientation preference and then analyzing its complexity."""

276 - def __call__(self,**params):

277 fm = super(measure_and_analyze_complexity,self).__call__(**params) 278 #from topo.command.analysis import measure_or_pref 279 #fm = measure_or_pref() 280 analyze_complexity(fm,simple_sheet_name="V1Simple",complex_sheet_name="V1Complex",filename="ModulationRatio")

281 282 283 pg= create_plotgroup(name='Orientation Preference and Complexity',category="Preference Maps", 284 doc='Measure preference for sine grating orientation.', 285 pre_plot_hooks=[measure_and_analyze_complexity.instance()]) 286 pg.add_plot('Orientation Preference',[('Hue','OrientationPreference')]) 287 pg.add_plot('Orientation Preference&Selectivity',[('Hue','OrientationPreference'), 288 ('Confidence','OrientationSelectivity')]) 289 pg.add_plot('Orientation Selectivity',[('Strength','OrientationSelectivity')]) 290 pg.add_plot('Modulation Ratio',[('Strength','ComplexSelectivity')]) 291 pg.add_plot('Phase Preference',[('Hue','PhasePreference')]) 292 pg.add_static_image('Color Key','command/or_key_white_vert_small.png') 293

Source Code for Module topo.analysis.vision