def triple_plot(cccsum, cccsum_hist, trace, threshold, save=False,
                savefile=''):
    r"""Main function to make a triple plot with a day-long seismogram, \
    day-long correlation sum trace and histogram of the correlation sum to \
    show normality.

    :type cccsum: numpy.ndarray
    :param cccsum: Array of the cross-channel cross-correlation sum
    :type cccsum_hist: numpy.ndarray
    :param cccsum_hist: cccsum for histogram plotting, can be the same as \
        cccsum but included if cccsum is just an envelope.
    :type trace: obspy.Trace
    :param trace: A sample trace from the same time as cccsum
    :type threshold: float
    :param threshold: Detection threshold within cccsum
    :type save: bool, optional
    :param save: If True will svae and not plot to screen, vice-versa if False
    :type savefile: str, optional
    :param savefile: Path to save figure to, only required if save=True
    """
    if len(cccsum) != len(trace.data):
        print('cccsum is: ' +
              str(len(cccsum))+' trace is: '+str(len(trace.data)))
        msg = ' '.join(['cccsum and trace must have the',
                        'same number of data points'])
        raise ValueError(msg)
    df = trace.stats.sampling_rate
    npts = trace.stats.npts
    t = np.arange(npts, dtype=np.float32) / (df * 3600)
    # Generate the subplot for the seismic data
    ax1 = plt.subplot2grid((2, 5), (0, 0), colspan=4)
    ax1.plot(t, trace.data, 'k')
    ax1.axis('tight')
    ax1.set_ylim([-15 * np.mean(np.abs(trace.data)),
                  15 * np.mean(np.abs(trace.data))])
    # Generate the subplot for the correlation sum data
    ax2 = plt.subplot2grid((2, 5), (1, 0), colspan=4, sharex=ax1)
    # Plot the threshold values
    ax2.plot([min(t), max(t)], [threshold, threshold], color='r', lw=1,
             label="Threshold")
    ax2.plot([min(t), max(t)], [-threshold, -threshold], color='r', lw=1)
    ax2.plot(t, cccsum, 'k')
    ax2.axis('tight')
    ax2.set_ylim([-1.7 * threshold, 1.7 * threshold])
    ax2.set_xlabel("Time after %s [hr]" % trace.stats.starttime.isoformat())
    # ax2.legend()
    # Generate a small subplot for the histogram of the cccsum data
    ax3 = plt.subplot2grid((2, 5), (1, 4), sharey=ax2)
    ax3.hist(cccsum_hist, 200, normed=1, histtype='stepfilled',
             orientation='horizontal', color='black')
    ax3.set_ylim([-5, 5])
    fig = plt.gcf()
    fig.suptitle(trace.id)
    fig.canvas.draw()
    if not save:
        plt.show()
        plt.close()
    else:
        plt.savefig(savefile)
    return

def survival_and_stats(feature, surv, upper_lim=5, axs=None, figsize=(7, 5), title=None,
                       order=None, colors=None, **args):
    if axs is None:
        fig = plt.figure(figsize=figsize)
        ax1 = plt.subplot2grid((3, 3), (0, 0), colspan=3, rowspan=2)
        ax2 = plt.subplot2grid((3, 3), (2, 0), colspan=2)
        ax3 = plt.subplot2grid((3, 3), (2, 2))
    else:
        ax1, ax2, ax3 = axs
        fig = plt.gcf()
    if feature.dtype != str:
        feature = feature.astype(str)
    if colors is None:
        colors = colors_global
    
    t = get_surv_fit(surv, feature)
    if order is None:
        t = t.sort([('5y Survival', 'Surv')], ascending=True)
    else:
        t = t.ix[order]
    survival_stat_plot(t, axs=[ax2, ax3], upper_lim=upper_lim, colors=colors)
    r = pd.Series({s:i for i, s in enumerate(t.index)})
    color_lookup = {c: colors[i % len(colors)] for i, c in enumerate(t.index)}
    
    draw_survival_curve(feature, surv, ax=ax1, colors=color_lookup, **args)
    ax1.legend().set_visible(False)
    if title:
        ax1.set_title(title)
    
    fig.tight_layout()

def plot_data(tag):
    data_array = tag.references[0]
    voltage = np.zeros(data_array.data.shape)
    data_array.data.read_direct(voltage)

    x_axis = data_array.dimensions[0]
    time = x_axis.axis(data_array.data_extent[0])

    spike_times = tag.positions[:]

    feature_data_array = tag.features[0].data
    snippets = tag.features[0].data[:]
    single_snippet = tag.retrieve_feature_data(3, 0)[:]

    snippet_time_dim = feature_data_array.dimensions[1]
    snippet_time = snippet_time_dim.axis(feature_data_array.data_extent[1])

    response_axis = plt.subplot2grid((2, 2), (0, 0), rowspan=1, colspan=2)
    single_snippet_axis = plt.subplot2grid((2, 2), (1, 0), rowspan=1, colspan=1)
    average_snippet_axis = plt.subplot2grid((2, 2), (1, 1), rowspan=1, colspan=1)

    response_axis.plot(time, voltage, color="dodgerblue", label=data_array.name)
    response_axis.scatter(spike_times, np.ones(spike_times.shape) * np.max(voltage), color="red", label=tag.name)
    response_axis.set_xlabel(x_axis.label + ((" [" + x_axis.unit + "]") if x_axis.unit else ""))
    response_axis.set_ylabel(data_array.label + ((" [" + data_array.unit + "]") if data_array.unit else ""))
    response_axis.set_title(data_array.name)
    response_axis.set_xlim(0, np.max(time))
    response_axis.set_ylim((1.2 * np.min(voltage), 1.2 * np.max(voltage)))
    response_axis.legend()

    single_snippet_axis.plot(snippet_time, single_snippet.T, color="red", label=("snippet No 4"))
    single_snippet_axis.set_xlabel(
        snippet_time_dim.label + ((" [" + snippet_time_dim.unit + "]") if snippet_time_dim.unit else "")
    )
    single_snippet_axis.set_ylabel(
        feature_data_array.label + ((" [" + feature_data_array.unit + "]") if feature_data_array.unit else "")
    )
    single_snippet_axis.set_title("single stimulus snippet")
    single_snippet_axis.set_xlim(np.min(snippet_time), np.max(snippet_time))
    single_snippet_axis.set_ylim((1.2 * np.min(snippets[3, :]), 1.2 * np.max(snippets[3, :])))
    single_snippet_axis.legend()

    mean_snippet = np.mean(snippets, axis=0)
    std_snippet = np.std(snippets, axis=0)
    average_snippet_axis.fill_between(
        snippet_time, mean_snippet + std_snippet, mean_snippet - std_snippet, color="red", alpha=0.5
    )
    average_snippet_axis.plot(snippet_time, mean_snippet, color="red", label=(feature_data_array.name + str(4)))
    average_snippet_axis.set_xlabel(
        snippet_time_dim.label + ((" [" + snippet_time_dim.unit + "]") if snippet_time_dim.unit else "")
    )
    average_snippet_axis.set_ylabel(
        feature_data_array.label + ((" [" + feature_data_array.unit + "]") if feature_data_array.unit else "")
    )
    average_snippet_axis.set_title("spike-triggered average")
    average_snippet_axis.set_xlim(np.min(snippet_time), np.max(snippet_time))
    average_snippet_axis.set_ylim((1.2 * np.min(mean_snippet - std_snippet), 1.2 * np.max(mean_snippet + std_snippet)))

    plt.subplots_adjust(left=0.15, top=0.875, bottom=0.1, right=0.98, hspace=0.35, wspace=0.25)
    plt.show()

def plot_state_timeseries(state, length):

    ax_main = pylab.subplot2grid((8, 8), (0, 0), colspan=4, rowspan=4)
    ax_main.grid(1)
    plot_state(ax_main, state, 10, length)

    for si, s in enumerate(["x", "y", "phi", "theta"]):
        ax_i = pylab.subplot2grid((8, 8), (4 + si, 0), colspan=8)
        ax_i.plot(state[s])
        ax_i.grid(1)

def createRegressionPlots(predictions,performance,coefs,fb_coefs,nfb_coefs,GroupDF,goodsubj,savefig=True):
    f=plt.figure(figsize=(22,12))
    ax1=plt.subplot2grid((2,4),(0,0), colspan=3)
    ax2=plt.subplot2grid((2,4),(0,3))
    ax3=plt.subplot2grid((2,4),(1,0), colspan=2)
    ax4=plt.subplot2grid((2,4),(1,2), colspan=2)

    dmnIdeal=pd.read_csv('/home/jmuraskin/Projects/NFB/analysis/DMN_ideal_2.csv')

    sns.tsplot(data=predictions,time='TR',value='predicted',unit='subj',condition='fb',ax=ax1)
    ax1.plot((dmnIdeal['Wander']-dmnIdeal['Focus'])/3,'k--')
    ax1.set_title('Average Predicted Time Series')

    g=sns.violinplot(data=performance,x='fb',y='R',split=True,bw=.3,inner='quartile',ax=ax2)
    # plt.close(g.fig)

    g=sns.violinplot(data=coefs,x='pe',y='Coef',hue='fb',split=True,bw=.3,inner='quartile',ax=ax3)
    g.plot([-1,len(unique(coefs['pe']))],[0,0],'k--')
    g.set_xlim([-.5,len(unique(coefs['pe']))])
    ylim=g.get_ylim()
    t,p = ttest_1samp(np.array(performance[performance.fb=='FEEDBACK']['R'])-np.array(performance[performance.fb=='NOFEEDBACK']['R']),0)
    ax2.set_title('Mean Subject Time Series Correlations-p=%0.2f' % p)

    t,p = ttest_1samp(np.array(fb_coefs['Coef'].reshape(len(unique(GroupDF[GroupDF.Subject_ID.isin(goodsubj)]['Subject_ID'])),len(unique(coefs['pe'])))),0)
    p05_FB,padj=fdr_correction(p,0.05)
    t,p = ttest_1samp(np.array(nfb_coefs['Coef'].reshape(len(unique(GroupDF[GroupDF.Subject_ID.isin(goodsubj)]['Subject_ID'])),len(unique(coefs['pe'])))),0)
    p05_NFB,padj=fdr_correction(p,0.05)
    for idx,(pFDR_FB,pFDR_NFB) in enumerate(zip(p05_FB,p05_NFB)):
        if pFDR_FB:
            ax3.scatter(idx,ylim[1]-.05,marker='*',s=75)
        if pFDR_NFB:
            ax3.scatter(idx,ylim[0]+.05,marker='*',s=75)


    t,p=ttest_1samp(np.array(fb_coefs['Coef']-nfb_coefs['Coef']).reshape(len(unique(GroupDF[GroupDF.Subject_ID.isin(goodsubj)]['Subject_ID'])),len(unique(coefs['pe']))),0)
    p05,padj=fdr_correction(p,0.05)

    sns.barplot(x=range(len(t)),y=t,ax=ax4,color='Red')
    for idx,pFDR in enumerate(p05):
        if pFDR:
            ax4.scatter(idx,t[idx]+ np.sign(t[idx])*0.2,marker='*',s=75)
    ax4.set_xlim([-0.5,len(unique(coefs['pe']))])
    ax4.set_xlabel('pe')
    ax4.set_ylabel('t-value')
    ax4.set_title('FB vs. nFB PE')

    for ax in [ax1,ax2,ax3,ax4]:
        for item in ([ax.title, ax.xaxis.label, ax.yaxis.label]):
            item.set_fontsize(18)
        for item in (ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(12)

    f.tight_layout()
    if savefig:
        f.savefig('%s/RSN_LinearRegPrediction.pdf' % saveFigureLocation,dpi=300)

def BS_plot(t, sol):
    def plot_p():
        ## Plot pump evolution
        fig_p, ax_p = pl.subplots() 
        for s in a_sol:
            ax_p.plot(t, s[:,2])
            ax_p.plot(t, s[:,3])

    ax_e = pl.subplot2grid((2,2),(1,0), colspan=2)
    ax_s = pl.subplot2grid((2,2),(0,0))
    ax_i = pl.subplot2grid((2,2),(0,1))
    
    plot_t(t, sol, (ax_s,ax_i))
    plot_e(t, sol, ax_e)

def plots1d(nimages=30, nreplicas=4, with_hist=True, show=True,
            height=5, width=8):
    plt.clf()
    fig = plt.gcf()
    fig.set_figheight(height)
    fig.set_figwidth(width)
    if with_hist:
        ncol_hist = 2
    else:
        ncol_hist = 0
    ax_list = [plt.subplot2grid((1,nimages+ncol_hist), (0,i), 
                               colspan=1, rowspan=nimages)
               for i in xrange(nimages)]
    
    for ax in ax_list:
        ax.set_ylim(0,1)
        ax.set_xticks([])
        ax.set_yticks([])
        ypos = np.random.uniform(0,1,nreplicas)
        ymax= ypos.max()
        ax.axhspan(0, ymax, alpha=0.2)
        xpos = np.zeros(nreplicas, )
        ax.scatter(xpos, ypos, c='k', facecolors="none")
        ax.scatter(0, ymax, c='r', linewidths=0, s=40)
    
    if with_hist:
        ax = plt.subplot2grid((1,nimages+ncol_hist), (0,nimages), 
                                   colspan=ncol_hist, rowspan=nimages)
        ax.set_xticks([])
        ax.set_yticks([])
        ax.set_ylim(0,1)

            
        n = 100000
        rmax = np.random.beta(nreplicas, 1, size=n)
        ax.hist(rmax, bins=np.sqrt(n)/10, orientation='horizontal', normed=True)
        
        if False:
            y = np.arange(1,0,-.01)
            x = y**(nreplicas-1)
            x /= x.max()
            ax.plot(x,y)
            ax.relim()
        
            
    
    
    if show:
        plt.show()

def plot_data(tag):
    data_array = tag.references[0]
    voltage = data_array[:]
    
    x_axis = data_array.dimensions[0]
    time = x_axis.axis(data_array.data_extent[0])
    
    stimulus_onset = tag.position
    stimulus_duration = tag.extent

    stimulus = tag.retrieve_feature_data(0)
    stimulus_array = tag.features[0].data
    
    stim_time_dim = stimulus_array.dimensions[0]
    stimulus_time = stim_time_dim.axis(stimulus_array.data_extent[0])
   
    response_axis = plt.subplot2grid((2, 2), (0, 0), rowspan=1, colspan=2)
    response_axis.tick_params(direction='out')
    response_axis.spines['top'].set_color('none')
    response_axis.spines['right'].set_color('none')
    response_axis.xaxis.set_ticks_position('bottom')
    response_axis.yaxis.set_ticks_position('left')

    stimulus_axis = plt.subplot2grid((2, 2), (1, 0), rowspan=1, colspan=2)
    stimulus_axis.tick_params(direction='out')
    stimulus_axis.spines['top'].set_color('none')
    stimulus_axis.spines['right'].set_color('none')
    stimulus_axis.xaxis.set_ticks_position('bottom')
    stimulus_axis.yaxis.set_ticks_position('left')
    
    response_axis.plot(time, voltage, color='dodgerblue', label=data_array.name, zorder=1)
    response_axis.set_xlabel(x_axis.label + ((" [" + x_axis.unit + "]") if x_axis.unit else ""))
    response_axis.set_ylabel(data_array.label + ((" [" + data_array.unit + "]") if data_array.unit else ""))
    response_axis.set_xlim(0, np.max(time))
    response_axis.set_ylim((1.2 * np.min(voltage), 1.2 * np.max(voltage)))
    response_axis.barh(1.2 * np.min(voltage), stim_duration, (1.2*np.max(voltage)) - (1.2*np.min(voltage)),
                       stim_onset, color='silver', alpha=0.25, zorder=0, label="stimulus epoch")
    response_axis.legend()
    
    stimulus_axis.plot(stimulus_time, stimulus, color="slategray", label="stimulus")
    stimulus_axis.set_xlabel(stim_time_dim.label + ((" [" + stim_time_dim.unit + "]") if stim_time_dim.unit else ""))
    stimulus_axis.set_ylabel(stimulus_array.label + ((" [" + stimulus_array.unit + "]") if stimulus_array.unit else ""))
    stimulus_axis.set_xlim(np.min(stimulus_time), np.max(stimulus_time))
    stimulus_axis.set_ylim(1.2 * np.min(stimulus), 1.2 * np.max(stimulus))
    stimulus_axis.legend()
    
    plt.subplots_adjust(left=0.15, top=0.875, bottom=0.1, right=0.98, hspace=0.45, wspace=0.25)
    plt.show()

def multiPlot(myData):    
    ax=plt.subplot2grid((ncol,nrow),(yPanel,xPanel))
    ax.hist(df[quantVar[myData]],alpha=0.5)#alpha??
    ax.set_xlabel(quantVar[myData],fontsize=14,fontweight='bold')
    #ax.set_ylabel('Y',fontsize=14,fontweight='bold',rotation='Vertical')
    plt.locator_params(axis = 'x', nbins = 2)
    ax.tick_params(labelsize=14)

def get_fig(height, n, layout='one_large_three_small'):
    
    import matplotlib.pylab as plt
    
    if layout == 'one_large_three_small':
    
        # Generate a figure
        fig = plt.figure(figsize = (4*height, 3*height))
        
        # Create the four subplots
        ax1 = plt.subplot2grid((3,4), (0,0), colspan=3, rowspan=3)
        ax2 = plt.subplot2grid((3,4), (0,3))
        ax3 = plt.subplot2grid((3,4), (1,3))
        ax4 = plt.subplot2grid((3,4), (2,3))

        # Create your axes list
        ax_list = [ ax1, ax2, ax3, ax4 ]

        fontsizes = [ 10*height, 5*height, 5*height, 5*height ] # Marker sizes
    

    if layout == 'one_row':
        # Generate a figure
        fig = plt.figure(figsize = (n*height, height))
        
        ax_list = []
        for i in range(n):
            ax = plt.subplot2grid((1,n), (0,i))
            ax_list.append(ax)
    
        fontsizes = [ 15, 15, 15, 15, 15, 15, 15, 15 ] # Marker sizes
         
         
    if layout == 'just_one':
        # Generate a figure
        fig = plt.figure(figsize = (height*1.5, height))

        # Add just one subplot
        ax1 = plt.axes()
        
        ax_list = [ ax1 ]
        
        fontsizes = [ 15 ]
            
    return fig, ax_list, fontsizes

def survival_stat_plot(t, upper_lim=5, axs=None, colors=None):
    """
    t is the DataFrame returned from a get_surv_fit call.
    """
    if axs is None:
        fig = plt.figure(figsize=(6, 1.5))
        ax = plt.subplot2grid((1, 3), (0, 0), colspan=2)
        ax2 = plt.subplot2grid((1, 3), (0, 2))
    else:
        ax, ax2 = axs
        fig = plt.gcf()
    if colors is None:
        colors = colors_global
    for i, (idx, v) in enumerate(t.iterrows()):
        conf_int = v['Median Survival']
        median_surv = v[('Median Survival', 'Median')]
        if (v['Stats']['# Events'] / v['Stats']['# Patients']) < .5:
            median_surv = np.nanmin([median_surv, 20])
            conf_int['Upper'] = np.nanmin([conf_int['Upper'], 20])
        l = ax.plot(*zip(*[[conf_int['Lower'], i], [median_surv, i], [conf_int['Upper'], i]]), lw=3, ls='--',
                    marker='o', dash_joinstyle='bevel', color=colors[i])
        ax.scatter(median_surv, i, marker='s', s=100, color=l[0].get_color(), edgecolors=['black'], zorder=10,
                   label=idx)
    ax.set_yticks(range(len(t)))
    ax.set_yticklabels(['{} ({})'.format(idx, int(t.ix[idx]['Stats']['# Patients'])) 
                        for idx in t.index])
    ax.set_ylim(-.5, i + .5)
    ax.set_xlim(0, upper_lim)
    ax.set_xlabel('Median Survival (Years)')
    
    tt = t['5y Survival']
    (tt['Surv']).plot(kind='barh', ax=ax2,
                      color=[l.get_color() for l in ax.lines],
                      xerr=[tt.Surv - tt.Lower, tt.Upper - tt.Surv],
                      width=.75,
                      ecolor='black')
    ax2.set_ybound(-.5,len(t)-.5)
    ax2.set_xlabel('5Y Survival')
    ax2.set_xticks([0, .5, 1.])
    ax2.set_yticks([])
    fig.tight_layout()

def plot_data(tag):
    data_array = tag.references[0]
    voltage = data_array[:]

    x_axis = data_array.dimensions[0]
    time = x_axis.axis(data_array.data_extent[0])

    spike_times = tag.positions[:]

    feature_data_array = tag.features[0].data
    stimulus = feature_data_array[:]

    stim_time_dim = feature_data_array.dimensions[0]
    stimulus_time = stim_time_dim.axis(feature_data_array.data_extent[0])

    response_axis = plt.subplot2grid((2, 2), (0, 0), rowspan=1, colspan=2)
    stimulus_axis = plt.subplot2grid((2, 2), (1, 0), rowspan=1, colspan=2, sharex=response_axis)

    response_axis.plot(time, voltage, color="dodgerblue", label=data_array.name)
    response_axis.scatter(spike_times, np.ones(spike_times.shape) * np.max(voltage), color="red", label=tag.name)
    response_axis.set_xlabel(x_axis.label + ((" [" + x_axis.unit + "]") if x_axis.unit else ""))
    response_axis.set_ylabel(data_array.label + ((" [" + data_array.unit + "]") if data_array.unit else ""))
    response_axis.set_title(data_array.name)
    response_axis.set_xlim(0, np.max(time))
    response_axis.set_ylim((1.2 * np.min(voltage), 1.2 * np.max(voltage)))
    response_axis.legend()

    stimulus_axis.plot(stimulus_time, stimulus, color="black", label="stimulus")
    stimulus_axis.scatter(spike_times, np.ones(spike_times.shape) * np.max(stimulus), color="red", label=tag.name)
    stimulus_axis.set_xlabel(stim_time_dim.label + ((" [" + stim_time_dim.unit + "]") if stim_time_dim.unit else ""))
    stimulus_axis.set_ylabel(
        feature_data_array.label + ((" [" + feature_data_array.unit + "]") if feature_data_array.unit else "")
    )
    stimulus_axis.set_title("stimulus")
    stimulus_axis.set_xlim(np.min(stimulus_time), np.max(stimulus_time))
    stimulus_axis.set_ylim(1.2 * np.min(stimulus), 1.2 * np.max(stimulus))
    stimulus_axis.legend()

    plt.subplots_adjust(left=0.15, top=0.875, bottom=0.1, right=0.98, hspace=0.45, wspace=0.25)
    # plt.savefig('taggedFeature.png')
    plt.show()

def get_fig(height, layout='one_large_three_small'):
    
    import matplotlib.pylab as plt
    
    if layout == 'one_large_three_small':
    
        # Generate a figure
        fig = plt.figure(figsize = (4*height, 3*height))
        
        # Create the four subplots
        ax1 = plt.subplot2grid((3,4), (0,0), colspan=3, rowspan=3)
        ax2 = plt.subplot2grid((3,4), (0,3))
        ax3 = plt.subplot2grid((3,4), (1,3))
        ax4 = plt.subplot2grid((3,4), (2,3))

        # Create your axes list
        ax = [ ax1, ax2, ax3, ax4 ]

        msizes = [ 50, 15, 15, 15 ] # Marker sizes
    
    if layout == 'four_equal_size_one_row':
        # Generate a figure
        fig = plt.figure(figsize = (4*height, height))
        
        # Create four subplots
        ax1 = plt.subplot2grid((1,4), (0,0))
        ax2 = plt.subplot2grid((1,4), (0,1))
        ax3 = plt.subplot2grid((1,4), (0,2))
        ax4 = plt.subplot2grid((1,4), (0,3))
    
        msizes = [ 5, 5, 5, 5 ] # Marker sizes
        
        # Create your axes list
        ax = [ ax1, ax2, ax3, ax4 ]
        
    if layout == 'just_one':
        # Generate a figure
        fig = plt.figure(figsize = (height*1.5, height))

        # Add just one subplot
        ax1 = plt.axes()
        
        ax = [ ax1 ]
        
        msizes = [ 5 ]
        
    # Set the marker sizes
    msizes = [ m * height for m in msizes ]
    
    return fig, ax, msizes

    def test_step(self):
        """Test function ``step``."""
        figure(); plot_shape = (1, 3)

        #Test SISO system
        A, B, C, D = self.make_SISO_mats()
        sys = ss(A, B, C, D)
        #print(sys)
        #print("gain:", dcgain(sys))

        subplot2grid(plot_shape, (0, 0))
        t, y = step(sys)
        plot(t, y)

        subplot2grid(plot_shape, (0, 1))
        T = linspace(0, 2, 100)
        X0 = array([1, 1])
        t, y = step(sys, T, X0)
        plot(t, y)

        #Test MIMO system
        A, B, C, D = self.make_MIMO_mats()
        sys = ss(A, B, C, D)

        subplot2grid(plot_shape, (0, 2))
        t, y = step(sys)
        plot(t, y)

def plot_update(fig, axarr, data1, data2, round):

    plt.cla()
    dota1 = data1 / STARTING_PCT
    dota2 = data2 / STARTING_PCT
    center_of_mass1 = ndimage.measurements.center_of_mass(dota1)
    center_of_mass2 = ndimage.measurements.center_of_mass(dota2)

    axarr = [plt.subplot(fig[0, 0]), plt.subplot2grid((2, 2), (1, 0), colspan=2), plt.subplot(fig[0, 1])]
    img1 = axarr[0].imshow(
        dota1,
        interpolation="nearest",
        cmap=plt.cm.ocean,
        extent=(0.5, np.shape(dota1)[0] + 0.5, 0.5, np.shape(dota1)[1] + 0.5),
    )
    axarr[0].plot([center_of_mass1[1]], [center_of_mass1[0]], "or")  # adds dot at center of mass
    axarr[0].plot([center_of_mass2[1]], [center_of_mass2[0]], "oy")  # adds dot for other teams c o m
    axarr[0].axis((1, 101, 1, 101))
    axarr[1].plot(avg_deal_data1, color="green")
    axarr[1].set_xlim(0, rounds)
    axarr[1].set_title("Current Round:" + str(round))
    axarr[0].set_title("Player1")
    axarr[2].set_title("Player2")
    axarr[0].set_ylabel("Give")
    axarr[0].set_xlabel("Accept")
    axarr[1].set_ylabel("Average Cash per Deal")
    axarr[1].set_xlabel("Round Number")
    plt.colorbar(img1, ax=axarr[0], label="Prevalence vs. Uniform")

    img2 = axarr[2].imshow(
        dota2,
        interpolation="nearest",
        cmap=plt.cm.ocean,
        extent=(0.5, np.shape(dota2)[0] + 0.5, 0.5, np.shape(dota2)[1] + 0.5),
    )
    axarr[2].plot([center_of_mass2[1]], [center_of_mass2[0]], "or")  # adds dot at center of mass
    axarr[2].plot([center_of_mass1[1]], [center_of_mass1[0]], "oy")  # adds dot for other teams c o m
    axarr[2].axis((1, 101, 1, 101))
    axarr[1].plot(avg_deal_data2, color="purple")

    plt.title("Current Round:" + str(round))
    axarr[2].set_ylabel("Give")
    axarr[2].set_xlabel("Accept")
    plt.colorbar(img2, ax=axarr[2], label="Prevalence vs. Uniform")

    plt.draw()

def plot_init(data1, data2):
    sns.set_style("dark")
    fig = gridspec.GridSpec(2, 2)
    axarr = [plt.subplot(fig[0, 0]), plt.subplot2grid((2, 2), (1, 0), colspan=2), plt.subplot(fig[0, 1])]
    img1 = axarr[0].imshow(
        data1,
        interpolation="nearest",
        cmap=plt.cm.ocean,
        extent=(0.5, np.shape(data1)[0] + 0.5, 0.5, np.shape(data1)[1] + 0.5),
    )
    plt.title("Current Round:" + str(0))
    axarr[0].set_ylabel("Give")
    axarr[0].set_xlabel("Accept")
    axarr[0].set_title("Distribution of Teams")
    axarr[1].plot(avg_deal_data1, color="green")
    axarr[1].set_xlim(0, rounds)
    axarr[1].set_ylabel("Average Cash per Deal")
    axarr[1].set_xlabel("Round Number")
    plt.colorbar(img1, ax=axarr[0], label="Prevalence vs. Uniform")

    img2 = axarr[2].imshow(
        data2,
        interpolation="nearest",
        cmap=plt.cm.ocean,
        extent=(0.5, np.shape(data2)[0] + 0.5, 0.5, np.shape(data2)[1] + 0.5),
    )
    plt.title("Current Round:" + str(0))
    axarr[2].set_ylabel("Give")
    axarr[2].set_xlabel("Accept")
    axarr[2].set_title("Distribution of Teams")
    axarr[1].plot(avg_deal_data2, color="purple")

    plt.colorbar(img2, ax=axarr[2], label="Prevalence vs. Uniform")

    plt.ion()
    mng = plt.get_current_fig_manager()
    mng.window.state("zoomed")
    plt.show()
    return fig, axarr

    def test_initial(self):
        A, B, C, D = self.make_SISO_mats()
        sys = ss(A, B, C, D)

        figure(); plot_shape = (1, 3)

        #X0=0 : must produce line at 0
        subplot2grid(plot_shape, (0, 0))
        t, y = initial(sys)
        plot(t, y)

        #X0=[1,1] : produces a spike
        subplot2grid(plot_shape, (0, 1))
        t, y = initial(sys, X0=matrix("1; 1"))
        plot(t, y)

        #Test MIMO system
        A, B, C, D = self.make_MIMO_mats()
        sys = ss(A, B, C, D)
        #X0=[1,1] : produces same spike as above spike
        subplot2grid(plot_shape, (0, 2))
        t, y = initial(sys, X0=[1, 1, 0, 0])
        plot(t, y)

#.........这里部分代码省略.........
    aperture_rad = math.ceil(float(fwhm_value)*2)      # Set aperture radius to three times the PSF radius
    sky_rad= math.ceil(aperture_rad)*5
    
    print aperture_rad, sky_rad

    if (not plot_only):

        if os.path.isfile(out_name): os.remove(out_name)
        if os.path.isfile(clean_name): os.remove(clean_name)

        # Check if files in list, otherwise exit
        if not ecfile:
           print "No .ec files in directory, exiting"
           sys.exit()
        
        
   
   
        iraf.noao.digiphot.apphot.qphot(image = image,\
        cbox = box ,\
        annulus = sky_rad ,\
        dannulus = 15. ,\
        aperture = str(aperture_rad),\
        coords = coords ,\
        output = out_name ,\
        plotfile = "" ,\
        zmag = 0. ,\
        exposure = "exptime" ,\
        airmass = "airmass" ,\
        filter = "filters" ,\
        obstime = "DATE" ,\
        epadu = gain ,\
        interactive = "no" ,\
        radplots = "yes" ,\
        verbose = "no" ,\
        graphics = "stdgraph" ,\
        display = "stdimage" ,\
        icommands = "" ,\
        wcsin = wcsin,
        wcsout = "logical",
        gcommands = "") 
        
         
        #iraf.noao.digiphot.apphot.phot(image=image, cbox=5., annulus=12.4, dannulus=10., salgori = "centroid", aperture=9.3,wcsin="world",wcsout="tv", interac = "no", coords=coords, output=out_name)
        iraf.txdump(out_name, "id,image,xcenter,ycenter,xshift,yshift,fwhm,msky,stdev,mag,merr", "yes", Stdout=clean_name)
        
    
    ma = np.genfromtxt(clean_name, comments="#", dtype=[("id","<f4"),  ("image","|S20"), ("X","<f4"), ("Y","<f4"), ("Xshift","<f4"), ("Yshift","<f4"),("fwhm","<f4"), ("ph_mag","<f4"), ("stdev","<f4"), ("fit_mag","<f4"), ("fiterr","<f4")])
    if (ma.size > 0):    
        m = ma[~np.isnan(ma["fit_mag"])]
    else:
        print "Only one object found!"
        m = np.array([ma])
        
    hdulist = pf.open(image)
    prihdr = hdulist[0].header
    img = hdulist[0].data * 1.
    nx, ny = img.shape

    
    
    dimX = int(4)
    dimY = int(np.ceil(len(m)*1./4))
    outerrad = sky_rad+10
    cutrad = outerrad + 15
    
    plt.suptitle("FWHM="+str(fwhm_value))
    k = 0
    for i in np.arange(dimX):
        for j in np.arange(dimY):
            if ( k < len(m)):
                ax = plt.subplot2grid((dimX,dimY),(i, j))
                y1, y2, x1, x2 = m[k]["X"]-cutrad, m[k]["X"]+cutrad, m[k]["Y"]-cutrad, m[k]["Y"]+cutrad
                y1, y2, x1, x2 = int(y1), int(y2), int(x1), int(x2)
                try:
                    zmin, zmax = zscale.zscale(img[x1:x2,y1:y2], nsamples=1000, contrast=0.25)
                except:
                    sh= img[x1:x2,y1:y2].shape
                    if sh[0]>0 and sh[1]>0:
                        zmin = np.nanmin(img[x1:x2,y1:y2])
                        zmax = np.nanmax(img[x1:x2,y1:y2])
                        continue
                    else:
                        continue
                ax.imshow(img[x1:x2,y1:y2], aspect="equal", extent=(-cutrad, cutrad, -cutrad, cutrad), origin="lower", cmap=matplotlib.cm.gray_r, interpolation="none", vmin=zmin, vmax=zmax)
                c1 = plt.Circle( (0, 0), edgecolor="r", facecolor="none", radius=5.)
                c2 = plt.Circle( (0, 0), edgecolor="orange", facecolor="none", radius=sky_rad)
                c3 = plt.Circle( (0, 0), edgecolor="yellow", facecolor="none", radius=sky_rad+10)
                plt.gca().add_artist(c1)
                plt.gca().add_artist(c2)
                plt.gca().add_artist(c3)
                ax.set_xticks([])
                ax.set_yticks([])
        
                plt.text(+5, +5, "%d"%m[k]["id"])
                plt.text(-cutrad, -cutrad, "%.2f$\pm$%.2f"%(m[k]["fit_mag"], m[k]["fiterr"]), color="b")
            k = k+1
    
    plt.savefig(os.path.join(plotdir, imname + "plot.png"))
    plt.clf()

def plot_pose(data, filename=None, gitsha1=None):
    t = pose.get_time_vector(data)
    ts = data['timestamp']
    dt = (ts - np.roll(ts, 1))[1:] # length is now 1 shorter than data

    # convert from system ticks to us
    t = t.astype('int') * 1000 * 1000 / 10000
    dt = dt.astype('int') * 1000 * 1000 / 10000

    rows = round((len(data.dtype.names) + 4) / 2)
    cols = 2
    names = data.dtype.names
    colors = sns.color_palette('husl', len(names))
    fig, axes = plt.subplots(rows, cols, sharex=True)
    axes = axes.ravel()

    base_title = 'bicycle pose'
    if filename is not None:
        base_title += ' (file \'{}\''.format(filename)
        if gitsha1 is not None:
            base_title += ', {}'.format(gitsha1)
        base_title += ')'

    title_size = mpl.rcParams['font.size'] + 2

    def set_title_func(f, x):
        title = base_title
        if x > 0:
            title = '{}\ndecimation factor {}'.format(base_title, x)
        f.suptitle(title, size=title_size)

    # convert angle data from radians to degrees
    angle_names = ('pitch', 'yaw', 'roll', 'steer', 'rear_wheel')
    data_np = []
    for name in names:
        if name in angle_names:
            data_np.append(data[name] * 180/np.pi)
        #elif name == 'timestamp':
        #    continue
        else:
            data_np.append(data[name])
        print('appending', name, 'to data')

    # plot objects to be added later
    dd_display = DecimatingDisplay(data_np, t, dt, set_title_func,
                                   None, None, None, None)

    full_time_range = (t.min(), t.max()) # start with full range of data
    td, datad, dfactor, indices = dd_display.decimate(full_time_range)
    set_title_func(fig, dfactor)

    lines = []
    for n, dd in enumerate(datad):
        axes_index = n + 4
        name = names[n]

        if name in angle_names:
            labelname = name + ' [°]'
        elif name in ('x', 'y'):
            labelname = name + ' [m]'
        elif name == 'v':
            labelname = name + ' [m/s]'
        elif name == 'timestamp':
            labelname = name + ' [system ticks @ 10 kHz]'
        elif name == 'computation_time':
            labelname = name + ' [us]'
        else:
            labelname = name
        print('creating line for', name)

        ax = axes[axes_index]
        lines.append(ax.plot(td, dd, label=labelname, color=colors[n])[0])
        ax.legend()
        ax.set_autoscale_on(False)
        ax.callbacks.connect('xlim_changed', dd_display.ax_update)
        if name == 'y':
            ax.invert_yaxis()

    # create dt plot
    ax = axes[3]
    lines.append(ax.plot(t[:-1], dt[indices[:-1]],
                 label='dt [us]', color=colors[-1])[0])
    ax.legend()
    ax.set_autoscale_on(False)
    ax.callbacks.connect('xlim_changed', dd_display.ax_update)

    axes[-1].set_xlabel('time [us]')
    axes[-2].set_xlabel('time [us]')

    # display trajectory plot
    ax = plt.subplot2grid((rows, cols), (0, 0), rowspan=2, sharey=axes[5])
    husl100 = mpl.colors.ListedColormap(sns.color_palette('husl', 100))
    lc, markers = _plot_trajectory(ax, data['x'][::100], data['y'][::100],
                                   t[::100], xpos='top', yinvert=True,
                                   cmap=husl100)
    axes[0] = ax # overwrite original axes

    # display histogram of sample time
    ax = plt.subplot2grid((rows, cols), (0, 1))
    histf = lambda dt: _plot_histogram(
#.........这里部分代码省略.........

 p.show()
 
 # Encode the filter output v=u*h with an IAF neuron
 # -------------------------------------------------
 # Specify parameters of the Ideal IAF neuron
 delta = 0.007                              # set the threshold
 bias  = 0.28                               # set the bias
 kappa = 1.                                 # set capacitance
 # # Encode the filter output
 spk_train, vol_trace = cim.iaf_encode(dt, t_proper, v_proper, b=bias, d=delta, C=kappa)
 run_time('Initialization time:',tic_init)  # display the initialization time 
 
 # Plot the voltage trace and the associated spike train
 p.figure(4)
 ax1 = p.subplot2grid((7, 1), (0, 0),rowspan=5,xlim=(0., 1.0),
         ylabel='Amplitude', ylim=(min(vol_trace), delta*1.1),
         title='Output of the [Filter]-[Ideal IAF] neural circuit')
 p.setp(ax1.get_xticklabels(), visible=False)
 ax2 = p.subplot2grid((7, 1), (5, 0),rowspan=2,sharex=ax1,
         xlabel='Time, [s]', ylim=(0., 1.1),xlim=(0., 1.0),
         yticks = (),yticklabels = ())
 ax1.plot(t_proper-t_proper[0],vol_trace,'b-',
          label='Membrane voltage $v\qquad$')
 ax1.plot((0, t_proper[-1]-t_proper[0]),(delta,delta),'--r',
          label='Threshold $\delta=' + repr(delta) + '$')
 ax1.plot(spk_train-t_proper[0],delta*np.ones_like(spk_train),'ro',
          label= '$v(t)=\delta$')
 ax1.legend(loc='lower right')
 ax2.stem(spk_train-t_proper[0],np.ones_like(spk_train),
          linefmt='k-', markerfmt='k^', label='$(t_k)_{k\in Z}$')
 ax2.legend(loc='lower right')