Plots and fixes, Smooth

- Fix: Show datapoint now uses coloured circles, just like the posterior plotting
- Added: Smooth function, taken from http://www.scipy.org/Cookbook/SignalSmooth
- Added: Pcolor_2d_data now refreshes data instead of replotting. Handles colorbar refresh as well. Fixed report_pixel accordingly
- Added: Function to do and plot the FFT power of some signal

datagenerator.py

@@ -10,6 +10,7 @@
 # from scaledimage import *
 import matplotlib.pyplot as plt
 import matplotlib.ticker as plttic
+import matplotlib.patches as plt_patches
 import numpy as np
 from scipy.spatial.distance import pdist
 
@@ -490,14 +491,12 @@ def show_datapoint_conjunctive(self, n=0):
         ax.set_yticklabels((r'$-\pi$', r'$-\frac{\pi}{2}$', r'$0$', r'$\frac{\pi}{2}$', r'$\pi$'))
 
         # Show ellipses at the stimuli positions
+        colmap = plt.get_cmap('gist_rainbow')
+        color_gen = [colmap(1.*(i)/self.T) for i in xrange(self.T)][::-1]  # use 22 colors
+
         for t in xrange(self.T):
-            e = Ellipse(xy=self.stimuli_correct[n, t], width=0.4, height=0.4)
-
-            ax.add_artist(e)
-            e.set_clip_box(ax.bbox)
-            e.set_alpha(0.5)
-            e.set_facecolor('white')
-            e.set_transform(ax.transData)
+            w = plt_patches.Wedge((self.stimuli_correct[n, t, 0], self.stimuli_correct[n, t, 1]), 0.25, 0, 360, 0.03, color=color_gen[t], alpha=0.7, linewidth=2)
+            ax.add_patch(w)
 
 
     def show_datapoint_features(self, n=0):

smooth.py

@@ -0,0 +1,104 @@
+##!/usr/bin/env python
+# encoding: utf-8
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def smooth(x, window_len=11, window='hanning'):
+    """smooth the data using a window with requested size.
+    
+    This method is based on the convolution of a scaled window with the signal.
+    The signal is prepared by introducing reflected copies of the signal 
+    (with the window size) in both ends so that transient parts are minimized
+    in the begining and end part of the output signal.
+    
+    input:
+        x: the input signal 
+        window_len: the dimension of the smoothing window; should be an odd integer
+        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
+            flat window will produce a moving average smoothing.
+
+    output:
+        the smoothed signal
+        
+    example:
+
+    t=linspace(-2,2,0.1)
+    x=sin(t)+randn(len(t))*0.1
+    y=smooth(x)
+    
+    see also: 
+    
+    np.hanning, np.hamming, np.bartlett, np.blackman, np.convolve
+    scipy.signal.lfilter
+ 
+    TODO: the window parameter could be the window itself if an array instead of a string
+    NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y.
+    """
+
+    if x.ndim != 1:
+        raise ValueError("smooth only accepts 1 dimension arrays.")
+
+    if x.size < window_len:
+        print x
+        raise ValueError("Input vector needs to be bigger than window size.")
+
+
+    if window_len<3:
+        return x
+
+
+    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
+        raise ValueError("Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'")
+
+
+    s=np.r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]
+    # print(len(s))
+    if window == 'flat':   # moving average
+        w=np.ones(window_len, 'd')
+    else:
+        w=eval('np.'+window+'(window_len)')
+
+    y=np.convolve(w/w.sum(), s, mode='valid')
+
+    # return y
+    return y[(window_len/2-1):-(window_len/2)]
+
+
+def smooth_demo():
+
+    t=np.linspace(-4, 4, 100)
+    x=np.sin(t)
+    xn=x+np.random.randn(len(t))*0.1
+    # y=smooth(x)
+
+    ws=31
+
+    plt.subplot(211)
+    plt.plot(np.ones(ws))
+
+    windows=['flat', 'hanning', 'hamming', 'bartlett', 'blackman']
+
+    plt.hold(True)
+    for w in windows[1:]:
+        eval('plt.plot(np.'+w+'(ws) )')
+
+    plt.axis([0, 30, 0, 1.1])
+
+    plt.legend(windows)
+    plt.title("The smoothing windows")
+    plt.subplot(212)
+    plt.plot(x)
+    plt.plot(xn)
+    for w in windows:
+        plt.plot(smooth(xn, 10, w))
+    l = ['original signal', 'signal with noise']
+    l.extend(windows)
+
+    plt.legend(l)
+    plt.title("Smoothing a noisy signal")
+    plt.show()
+
+
+if __name__=='__main__':
+    smooth_demo()

utils.py

@@ -26,7 +26,7 @@
 import datetime
 
 __maxexp__ = np.finfo('float').maxexp
-
+plt.ion()
 
 ############################ MATH ##################################
 
@@ -585,71 +585,111 @@ def pcolor_2d_data(data, x=None, y=None, xlabel='', ylabel='', title='', colorba
         f = plt.figure(fignum)
         ax_handle = f.add_subplot(111)
         ax_handle.clear()
-
-    if log_scale:
-        im = ax_handle.imshow(data.T, interpolation=interpolation, origin='lower left', norm=LogNorm())
     else:
-        im = ax_handle.imshow(data.T, interpolation=interpolation, origin='lower left')
+        plt.figure(ax_handle.get_figure().number)
 
-    if not x is None:
-        assert data.shape[0] == x.size, 'Wrong x dimension'
+    if len(ax_handle.get_images()) > 0:
+        # Update the data if the figure is already filled
 
-        if not ticks_interpolate is None:
-            selected_ticks = np.array(np.linspace(0, x.size-1, ticks_interpolate), dtype=int)
-            ax_handle.set_xticks(selected_ticks)
-            ax_handle.set_xticklabels([label_format % x[tick_i] for tick_i in selected_ticks], rotation=90)
-        else:
-            ax_handle.set_xticks(np.arange(x.size))
-            ax_handle.set_xticklabels([label_format % curr for curr in x], rotation=90)
+        im = ax_handle.get_images()[0]
+        im.set_data(data.T)
+        im.set_clim(vmin=np.nanmin(data), vmax=np.nanmax(data))
+        im.changed()
 
-    if not y is None:
-        assert data.shape[1] == y.size, 'Wrong y dimension'
+        # Change mouse over behaviour
+        def report_pixel(x_mouse, y_mouse): 
+            # Extract loglik at that position
 
-        if not ticks_interpolate is None:
-            selected_ticks = np.array(np.linspace(0, y.size-1, ticks_interpolate), dtype=int)
-            ax_handle.set_yticks(selected_ticks)
-            ax_handle.set_yticklabels([label_format % y[tick_i] for tick_i in selected_ticks])
-        else:
-            ax_handle.set_yticks(np.arange(y.size))
-            ax_handle.set_yticklabels([label_format % curr for curr in y])
-
-    if xlabel:
-        ax_handle.set_xlabel(xlabel)
+            try:
+                x_i = int(np.round(x_mouse))
+                y_i = int(np.round(y_mouse))
+
+                if x is not None:
+                    x_display = x[x_i]
+                else:
+                    x_display = x_i
+
+                if y is not None:
+                    y_display = y[y_i]
+                else:
+                    y_display = y_i
+
+                return "x=%.2f y=%.2f value=%.2f" % (x_display, y_display, data[x_i, y_i])
+            except:
+                return ""
+
+        ax_handle.format_coord = report_pixel
 
-    if ylabel:
-        ax_handle.set_ylabel(ylabel)
 
-    if colorbar:
-        ax_handle.get_figure().colorbar(im)
-
-    if title:
-        ax_handle.set_title(title)
-
-    ax_handle.axis('tight')
+    else:
+        # Create the Figure
+        if log_scale:
+            im = ax_handle.imshow(data.T, interpolation=interpolation, origin='lower left', norm=LogNorm())
+        else:
+            im = ax_handle.imshow(data.T, interpolation=interpolation, origin='lower left')
 
-    # Change mouse over behaviour
-    def report_pixel(x_mouse, y_mouse): 
-        # Extract loglik at that position
+        if not x is None:
+            assert data.shape[0] == x.size, 'Wrong x dimension'
 
-        try:
-            x_i = int(x_mouse)
-            y_i = int(y_mouse)
-
-            if x is not None:
-                x_display = x[x_i]
+            if not ticks_interpolate is None:
+                selected_ticks = np.array(np.linspace(0, x.size-1, ticks_interpolate), dtype=int)
+                ax_handle.set_xticks(selected_ticks)
+                ax_handle.set_xticklabels([label_format % x[tick_i] for tick_i in selected_ticks], rotation=90)
             else:
-                x_display = x_i
-
-            if y is not None:
-                y_display = y[y_i]
+                ax_handle.set_xticks(np.arange(x.size))
+                ax_handle.set_xticklabels([label_format % curr for curr in x], rotation=90)
+
+        if not y is None:
+            assert data.shape[1] == y.size, 'Wrong y dimension'
+
+            if not ticks_interpolate is None:
+                selected_ticks = np.array(np.linspace(0, y.size-1, ticks_interpolate), dtype=int)
+                ax_handle.set_yticks(selected_ticks)
+                ax_handle.set_yticklabels([label_format % y[tick_i] for tick_i in selected_ticks])
             else:
-                y_display = y_i
+                ax_handle.set_yticks(np.arange(y.size))
+                ax_handle.set_yticklabels([label_format % curr for curr in y])
 
-            return "x=%.2f y=%.2f value=%.2f" % (x_display, y_display, data[x_i, y_i])
-        except:
-            return ""
+        if xlabel:
+            ax_handle.set_xlabel(xlabel)
+
+        if ylabel:
+            ax_handle.set_ylabel(ylabel)
+
+        if colorbar:
+            ax_handle.get_figure().colorbar(im)
+
+        if title:
+            ax_handle.set_title(title)
 
-    ax_handle.format_coord = report_pixel
+        ax_handle.axis('tight')
+
+        # Change mouse over behaviour
+        def report_pixel(x_mouse, y_mouse): 
+            # Extract loglik at that position
+
+            try:
+                x_i = int(np.round(x_mouse))
+                y_i = int(np.round(y_mouse))
+
+                if x is not None:
+                    x_display = x[x_i]
+                else:
+                    x_display = x_i
+
+                if y is not None:
+                    y_display = y[y_i]
+                else:
+                    y_display = y_i
+
+                return "x=%.2f y=%.2f value=%.2f" % (x_display, y_display, data[x_i, y_i])
+            except:
+                return ""
+
+        ax_handle.format_coord = report_pixel
+
+    # redraw
+    ax_handle.get_figure().canvas.draw()
 
     return ax_handle, im
 
@@ -800,7 +840,7 @@ def plot_torus(theta, gamma, Z, weight_deform=0., torus_radius=5., tube_radius=3
             cb = mplt.colorbar(title='', orientation='vertical', label_fmt='%.2f', nb_labels=5)
 
         mplt.outline(color=(0., 0., 0.))
-        mplt.show()
+        mplt.draw()
 
     else:
         fig = plt.figure()
@@ -815,7 +855,7 @@ def plot_torus(theta, gamma, Z, weight_deform=0., torus_radius=5., tube_radius=3
         if draw_colorbar:
             plt.colorbar(m)
 
-        plt.show()
+        # plt.show()
 
 
 def plot_powerlaw_fit(xdata, ydata, amp, index, yerr=None, fignum=None):
@@ -850,6 +890,20 @@ def plot_powerlaw_fit(xdata, ydata, amp, index, yerr=None, fignum=None):
     plt.xlim((xdata.min()*0.9, xdata.max()*1.1))
 
 
+def plot_fft_power(data, dt=1.0, n=None, axis=-1, fignum=None):
+    '''
+        Compute the FFT and plot the power spectrum.
+    '''
+
+    freq = np.fft.fftfreq(data.shape[-1], d=dt)
+    FS = np.fft.fft(data, n=n, axis=axis)
+
+    plt.figure(fignum)
+    # plt.plot(freq, np.log(np.abs(np.fft.fftshift(FS))**2.))
+    plt.plot(freq, np.abs(FS)**2.)
+    plt.title('Power spectrum')
+
+
 def plot_fft2_power(data, s=None, axes=(-2, -1), fignum=None):
     '''
         Compute the 2D FFT and plot the 2D power spectrum.