def plot_x_query(lwr_theta):
    print 'Welcome to the plot demo!\n================================\n Locally weighted regression: lwr_theta for an x_query:\n%s' % lwr_theta

    # Read data to make a scatter plot and LR fit
    x, y = read_data('/nosql/input/x_data.txt')
    x_2dim = x[:, [0, COLUMN_NUMBER]]  # because we can't plot more than two dimensions
    theta = batch_linear_regression(x_2dim, y)
    print 'Simple liner regression produce the following theta for the regression line:\n%s' % theta

    # construct a line from intercept and coefficient
    x_pred = [min(x[:, 1])[0, 0], max(x[:, 1])[0, 0]]
    y_pred = [theta[0, 0] + v * theta[1, 0] for v in x_pred]

    plt.figure(figsize=(14,10))
    plt.plot(x_2dim[:, 1], y, 'bo', label='Training Set')
    plt.plot(x_pred, y_pred, 'r-', label='Linear Regression')

    # Predict outcome for x_query
    y_query = x_query * lwr_theta
    print 'Given the x_query: %s\nLWR predicts target: %s' % (x_query, y_query)
    plt.plot(x_query[:, COLUMN_NUMBER], y_query, 'go', markersize=10,
            label='Locally Wheighted Linear Regression x_query Prediction')

    # Circle the prediction and fine tune the plot
    circle = plt.Circle((x_query[:, COLUMN_NUMBER], y_query), 2, color='y', fill=False)
    plt.gca().add_artist(circle)
    plt.grid()
    plt.legend(loc=2)
    plt.tight_layout()
    plt.show()

def map_along_line(x, y, q, ax=None, cmap=None, norm=None,
            time=None, max_step=1., missing=np.nan,
            new_timebase=None,
            **kwargs):
    """Map some quantity q along x,y as a coloured line.
With time set, perform linear interpolation of x,y,q onto new_timebase
filling with missing, and with max_step."""

    if ax is None:
        ax = plt.gca()

    if x.shape != y.shape:
        raise ValueError('Shape mismatch')
    if x.shape != q.shape:
        raise ValueError('Shape mismatch')

    if time is not None:
        if new_timebase is None:
            new_timebase = np.arange(time[0], time[-1], np.min(np.diff(time)))

        # Bit redundant
        x = interp_safe(new_timebase, time, x, max_step=max_step, missing=missing)
        y = interp_safe(new_timebase, time, y, max_step=max_step, missing=missing)
        q = interp_safe(new_timebase, time, q, max_step=max_step, missing=missing)

    points = np.array([x, y]).T.reshape(-1, 1, 2)
    segments = np.concatenate([points[:-1], points[1:]], axis=1)
    lc = LineCollection(segments, cmap=cmap, norm=norm, **kwargs)

    lc.set_array(q)
    plt.gca().add_collection(lc)

    return lc

def plot_svc(X, y, mysvc, bounds=None, grid=50):
    if bounds is None:
        xmin = np.min(X[:, 0], 0)
        xmax = np.max(X[:, 0], 0)
        ymin = np.min(X[:, 1], 0)
        ymax = np.max(X[:, 1], 0)
    else:
        xmin, ymin = bounds[0], bounds[0]
        xmax, ymax = bounds[1], bounds[1]
    aspect_ratio = (xmax - xmin) / (ymax - ymin)
    xgrid, ygrid = np.meshgrid(np.linspace(xmin, xmax, grid),
                              np.linspace(ymin, ymax, grid))
    plt.gca(aspect=aspect_ratio)
    plt.xlim(xmin, xmax)
    plt.ylim(ymin, ymax)
    plt.xticks([])
    plt.yticks([])
    plt.hold(True)
    plt.plot(X[y == 1, 0], X[y == 1, 1], 'bo')
    plt.plot(X[y == -1, 0], X[y == -1, 1], 'ro')
    
    box_xy = np.append(xgrid.reshape(xgrid.size, 1), ygrid.reshape(ygrid.size, 1), 1)
    if mysvc is not None:
        scores = mysvc.decision_function(box_xy)
    else:
        print 'You must have a valid SVC object.'
        return None;
    
    CS=plt.contourf(xgrid, ygrid, scores.reshape(xgrid.shape), alpha=0.5, cmap='jet_r')
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[0], colors='k', linestyles='solid', linewidths=1.5)
    plt.contour(xgrid, ygrid, scores.reshape(xgrid.shape), levels=[-1,1], colors='k', linestyles='dashed', linewidths=1)
    plt.plot(mysvc.support_vectors_[:,0], mysvc.support_vectors_[:,1], 'ko', markerfacecolor='none', markersize=10)
    CB = plt.colorbar(CS)

    def on_circle_click(self, x, y):
        if not self.scale["length_of_scale_in_mu_meter"]:
            self.on_scala_click(x, y)
            return
        plt.scatter(x, y, color="white", s=60, edgecolor="blue", marker="x",
            lw=2.5)

        if len(self.current_circle_points) == 3:
            self.current_circle_points[:] = []
        self.current_circle_points.append((x, y))
        if len(self.current_circle_points) == 3:
            p = self.current_circle_points

            point, radius = cirle_from_three_points(p[0][0], p[0][1], p[1][0],
                p[1][1], p[2][0], p[2][1])

            color = self.colors[self.color_index % len(self.colors)]
            self.color_index += 1

            circle = matplotlib.patches.Circle(point, radius, lw=4,
                facecolor="none", edgecolor=color)

            factor = self.scale["length_of_scale_in_mu_meter"] / \
                self.scale["length_of_scale_in_px"]

            plt.text(point[0], point[1], "Radius: %.4fE-6 m" % (radius *
                factor), horizontalalignment='center', color="black",
                verticalalignment='center', fontsize=11,
                bbox=dict(facecolor=color, alpha=0.85, edgecolor="0.7"))

            plt.gca().add_patch(circle)
        self.redraw()

def mwpc_last_spills(df, filename, n_spills):
    spacing = 6.0 # mm
    L = spacing*len(df.T)
    y_max = df.max().max()

    fig = plt.figure()
    ax = fig.add_subplot(111)
    x = ((np.array(df.sum().index))*spacing)-(L/2.0)+(spacing/2)

    colormap = plt.cm.gist_ncar
    plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, n_spills)])

    for i, index in enumerate(df.index):
        y = df.ix[i]
        index = str(index)
        stime = index.split(' ')[1].split('.')[0]
        sdate = index.split(' ')[0]
        ax.plot(x,y, ls='steps-mid', label=stime, linewidth=2)

    ax.grid()
    ax.set_xlabel('x [mm]')
    ax.set_ylabel('y')
    ax.set_title('MWPC ({}) {}'.format(sdate, filename))
    ax.legend(fancybox=True, framealpha=0.5)
    ax.set_ylim(top=y_max*1.2)

    figname = filename.split('_')[1].split('.')[0]
    filename_figure = '{}{}_1.png'.format(directory, figname)
    try:
        fig.savefig(filename_figure)
    except IOError:
        print('{} access denied!'.format(filename_figure))
        pass
    plt.close(fig)
    return df

def mass_flux_plot(*args,**kwargs):
    fltm = idls.read(args[0])
    injm = idls.read(args[1])
    f1 = plt.figure()

    ax1 = f1.add_subplot(211)
    plt.plot(injm.nt_sc,injm.nmf_rscale,'r')
    plt.plot(injm.nt_sc,injm.nmf_zscale,'b')
    plt.plot(injm.nt_sc,injm.nmf_z0scale,'k')
    plt.plot(injm.nt_sc,(injm.nmf_rscale+injm.nmf_zscale),'g')
    plt.axis([0.0,160.0,0.0,3.5e-5])
    plt.minorticks_on()
    locs,labels = plt.yticks()
    plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e5))
    plt.text(0.0, 1.03, r'$10^{-5}$', transform = plt.gca().transAxes)
    plt.xlabel(r'Time [yr]',labelpad=6)
    plt.ylabel(r'$\dot{\rm M}_{\rm out} [ \rm{M}_{\odot} \rm{yr}^{-1} ]$',labelpad=15)
    
    ax2 = f1.add_subplot(212)
    plt.plot(fltm.nt_sc,fltm.nmf_rscale,'r')
    plt.plot(fltm.nt_sc,fltm.nmf_zscale,'b')
    plt.plot(fltm.nt_sc,fltm.nmf_z0scale,'k')
    plt.plot(fltm.nt_sc,(fltm.nmf_rscale+fltm.nmf_zscale),'g')
    plt.axis([0.0,160.0,0.0,4.0e-5])
    plt.minorticks_on()
    locs,labels = plt.yticks()
    plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e5))
    plt.text(0.0, 1.03, r'$10^{-5}$', transform = plt.gca().transAxes)
    plt.xlabel(r'Time [yr]',labelpad=6)
    plt.ylabel(r'$\dot{\rm M}_{\rm out} [ \rm {M}_{\odot} \rm{yr}^{-1} ]$',labelpad=15)

def graphical_test(satisfactory=0):
    from matplotlib import cm, pylab
    def cons():
        return np.random.random(2)*4-2

    def foo(x,y,a,b):
        "banana function"
        tmp=a-x
        tmp*=tmp
        out=-x*x
        out+=y
        out*=out
        out*=b
        out+=tmp
        return out*(abs(np.cos((x-1)**2+(y-1)**2))+10.0/b)
    def f(params):
        return foo(params[0], params[1],1,100)

    optimizer=optimize(f, cons, verbose=False,its=1, hillWalks=0, satisfactory=satisfactory, finalWalk=0)
    
    bgx,bgy=np.mgrid[-2:2:1000j,-2:2:1000j]
    bg=foo(bgx,bgy, 1,100)
    for i in xrange(20):
        pylab.clf()
        pylab.imshow(bg, cmap=cm.RdBu,vmax=bg.mean()/10)
        for x in optimizer.pool: pylab.plot((x[2]+2)/4*1000,(x[1]+2)/4*1000, ('gx'))
        print optimizer.pool[0],optimizer.muterate
        pylab.gca().set_xbound(0,1000)
        pylab.gca().set_ybound(0,1000)
        pylab.draw()
        pylab.colorbar()
        optimizer.run()
        raw_input('enter to advance')
    return optimizer

def AR_predict(d,coeffs, names=None, plot=False):
	"""
	Plot the auto-regression predictions and the actual data.
	"""
	predictions, ax1, ax2 = [], None, None
	i = 0
	for c in coeffs:
		p = len(c)
		y_predict = np.convolve(d,c[::-1],mode='valid')
		y_predict = y_predict[:-1] ## discard the last value because its outside our domain
		predictions.append(y_predict)

		if plot:
			series_name = names[i] if names!=None else ""
			y_gt = d[p:]
			N = len(y_gt)
			plt.subplot(2,1,1)
			if ax1== None:
				ax1 = plt.gca()
				ax1.plot(np.arange(N), y_gt, label="actual")
			ax1.plot(np.arange(N), y_predict, label="prediction %s (p=%d)"%(series_name,p))
			ax1.legend()
			plt.subplot(2,1,2)
			if ax2==None: ax2 = plt.gca()
			ax2.plot(np.arange(p), c[::-1], label= series_name+' coefficients')
			ax2.legend()
			i += 1
	if plot: plt.show()

	return predictions

def _gaussian_test():
    import matplotlib.pyplot as plt
    n = 10000
    mu_x = 0.0
    mu_y = 0.0
    #sig_x, sig_y = 1.5, 1.5
    tau = 0.0
    seeing = 1.5
    sigma = seeing / (2. * np.sqrt(2. * np.e))
    slit_width = 0.2
    slit_height = 10.0
    slit_x = np.empty(n, dtype=np.float64)
    slit_y = np.empty(n, dtype=np.float64)
    slit_x, slit_y = slit_gaussian_psf(n, mu_x, mu_y, sigma, sigma, tau, slit_width, slit_height)
    log.info("x range: [%s, %s]", slit_x.min(), slit_x.max())
    log.info("y range: [%s, %s]", slit_y.min(), slit_y.max())
    plt.scatter(slit_x, slit_y, alpha=0.8)
    plt.fill([-slit_width/2, slit_width/2, slit_width/2, -slit_width/2],
             [-slit_height/2, -slit_height/2, slit_height/2, slit_height/2],
             'r',
             alpha=0.10,
             edgecolor='k')
    plt.gca().set_aspect("equal")
    plt.title("Gaussian distribution")
    plt.xlim([-slit_height/2., slit_height/2])
    plt.show()

    def check_models(self):
        temp = np.logspace(0, np.log10(600))
        num = len(self.available_models())

        fig, ax = plt.subplots(1)
        self.plotting_colours(num, fig, ax, repeats=2)

        for author in self.available_models():
            Nc, Nv = self.update(temp=temp, author=author)
            # print Nc.shape, Nv.shape, temp.shape
            ax.plot(temp, Nc, '--')
            ax.plot(temp, Nv, '.', label=author)

        ax.loglog()
        leg1 = ax.legend(loc=0, title='colour legend')

        Nc, = ax.plot(np.inf, np.inf, 'k--', label='Nc')
        Nv, = ax.plot(np.inf, np.inf, 'k.', label='Nv')

        plt.legend([Nc, Nv], ['Nc', 'Nv'], loc=4, title='Line legend')
        plt.gca().add_artist(leg1)

        ax.set_xlabel('Temperature (K)')
        ax.set_ylabel('Density of states (cm$^{-3}$)')
        plt.show()

    def plot_timeseries(self, ax=None, vmin=None, vmax=None,
            colorbar=False, label=True):

        if vmin is None:
            vmin = self.vmin

        if vmax is None:
            vmax = self.vmax

        if ax is None:
            ax = plt.gca()
        plt.sca(ax)
        plt.cla()
        plt.imshow(self.tser_arr[::-1,:], vmin=vmin, vmax=vmax,
            interpolation='Nearest', extent=self.extent, origin='upper',aspect='auto')
        plt.xlim(self.extent[0], self.extent[1])
        plt.ylim(self.extent[2], self.extent[3])
        # plt.vlines(self.ionogram_list[0].time, self.extent[2], self.extent[3], 'r')
        if label:
           celsius.ylabel('f / MHz')

        if colorbar:
            old_ax = plt.gca()
            plt.colorbar(
                    cax = celsius.make_colorbar_cax(), ticks=self.cbar_ticks
                ).set_label(r"$Log_{10} V^2 m^{-2} Hz^{-1}$")
            plt.sca(old_ax)

    def __init__(self, history_length=100):

        self.history_length = history_length

        self.root = Tk.Tk()
        self.root.wm_title("GPU MEMORY DISPLAY")
        self.root.protocol('WM_DELETE_WINDOW', self.quit_button_cb)

        fig = plt.figure(figsize=(8, 3))

        self.subplot = plt.subplot(211)
        plt.gca().invert_xaxis()

        self.canvas = FigureCanvasTkAgg(fig, master=self.root)

        self.max_gpu_mem = None
        self.current_gpu_mem = None
        self.mem_data = [0] * self.history_length
        self.mem_range = list(reversed(range(self.history_length)))

        self.canvas.get_tk_widget().pack(side=Tk.TOP, fill=Tk.BOTH, expand=1)

        self.update()

        Tk.mainloop()

    def showladderpca(self):

        import mdp
        from matplotlib import pylab as plt
        import pprint
        import math

        cerr('calculating PCA & plotting')
        peak_sizes = sorted(list([ x.rtime for x in self.alleles ]))
        #peak_sizes = sorted( peak_sizes )[:-5]
        #pprint.pprint(peak_sizes)
        #comps = algo.simple_pca( peak_sizes )
        #algo.plot_pca(comps, peak_sizes)

        from fatools.lib import const
        std_sizes = const.ladders['LIZ600']['sizes']

        x = std_sizes
        y = [ x * 0.1 for x in peak_sizes ]

        D = np.zeros( (len(y), len(x)) )
        for i in range(len(y)):
            for j in range(len(x)):
                D[i,j] = math.exp( ((x[j] - y[i]) * 0.001) ** 2 )

        pprint.pprint(D)
        im = plt.imshow(D, interpolation='nearest', cmap='Reds')
        plt.gca().invert_yaxis()
        plt.xlabel("STD")
        plt.ylabel("PEAK")
        plt.grid()
        plt.colorbar()
        plt.show()

    def plot_std_meshlines(self, step=0.1):
        '''
        plot mesh circles for stdv
        '''

        color = self.std_color

        nstdmax = self.stdmax
        if self.negative:
            axmin = -np.pi / 2.
        else:
            axmin = 0.

        th = np.arange(axmin, np.pi / 2, 0.01)

        for ra in np.arange(0, nstdmax + 0.1 * step, step):
            self.ax.plot(ra * np.sin(th), ra * np.cos(th), ':', color=color)

        if self.normalize:
            self.ax.set_ylabel('$\sigma / \sigma_{obs}$', color=color)
            self.ax.set_xlabel('$\sigma / \sigma_{obs}$', color=color)
        else:
            self.ax.set_ylabel('Standard Deviation', color=color)
            self.ax.set_xlabel('Standard Deviation', color=color)

        xticklabels = plt.getp(plt.gca(), 'xticklabels')
        plt.setp(xticklabels, color=color)
        yticklabels = plt.getp(plt.gca(), 'yticklabels')
        plt.setp(yticklabels, color=color)

def hideaxis(pos=None):
    # hide x y axis
    if pos:
        df = pd.DataFrame(pos.values(), columns=['x', 'y'])
        plt.xlim([df['x'].min()-5, df['x'].max()+5])
        plt.ylim([df['y'].min()-5, df['y'].max()+5])
    plt.gca().xaxis.set_major_locator(plt.NullLocator())
    plt.gca().yaxis.set_major_locator(plt.NullLocator())

def dateticks(fmt='%Y-%m', **kwargs):
    '''setup the date ticks'''
    dateticker = ticker.FuncFormatter(lambda numdate, _: num2date(numdate).strftime(fmt))
    pylab.gca().xaxis.set_major_formatter(dateticker)
    # pylab.gcf().autofmt_xdate()
    tmp = dict(rotation=30, ha='right')
    tmp.update(kwargs)
    pylab.setp(pylab.xticks()[1], **tmp)

def continuous_calibration():
    utc = Calendar()
    t_start = utc.time(YMDhms(2011, 9, 1))
    t_fc_start = utc.time(YMDhms(2015, 10, 1))
    dt = deltahours(1)
    n_obs = int(round((t_fc_start - t_start)/dt))
    obs_time_axis = Timeaxis(t_start, dt, n_obs + 1)
    q_obs_m3s_ts = observed_tistel_discharge(obs_time_axis.total_period())

    ptgsk = create_tistel_simulator(PTGSKOptModel, tistel.geo_ts_repository(tistel.grid_spec.epsg()))
    initial_state = burn_in_state(ptgsk, t_start, utc.time(YMDhms(2012, 9, 1)), q_obs_m3s_ts)

    num_opt_days = 30
    # Step forward num_opt_days days and store the state for each day:
    recal_start = t_start + deltahours(num_opt_days*24)
    t = t_start
    state = initial_state
    opt_states = {t: state}
    while t < recal_start:
        ptgsk.run(Timeaxis(t, dt, 24), state)
        t += deltahours(24)
        state = ptgsk.reg_model_state
        opt_states[t] = state

    recal_stop = utc.time(YMDhms(2011, 10, 30))
    recal_stop = utc.time(YMDhms(2012, 5, 30))
    curr_time = recal_start
    q_obs_avg = TsTransform().to_average(t_start, dt, n_obs + 1, q_obs_m3s_ts)
    target_spec = TargetSpecificationPts(q_obs_avg, IntVector([0]), 1.0, KLING_GUPTA)
    target_spec_vec = TargetSpecificationVector([target_spec])
    i = 0
    times = []
    values = []
    p, p_min, p_max = construct_calibration_parameters(ptgsk)
    while curr_time < recal_stop:
        print(i)
        i += 1
        opt_start = curr_time - deltahours(24*num_opt_days)
        opt_state = opt_states.pop(opt_start)
        p = ptgsk.region_model.get_region_parameter()
        p_opt = ptgsk.optimize(Timeaxis(opt_start, dt, 24*num_opt_days), opt_state, target_spec_vec,
                               p, p_min, p_max, tr_stop=1.0e-5)
        ptgsk.region_model.set_region_parameter(p_opt)
        corr_state = adjust_simulator_state(ptgsk, curr_time, q_obs_m3s_ts)
        ptgsk.run(Timeaxis(curr_time, dt, 24), corr_state)
        curr_time += deltahours(24)
        opt_states[curr_time] = ptgsk.reg_model_state
        discharge = ptgsk.region_model.statistics.discharge([0])
        times.extend(discharge.time(i) for i in range(discharge.size()))
        values.extend(list(np.array(discharge.v)))
    plt.plot(utc_to_greg(times), values)
    plot_results(None, q_obs=observed_tistel_discharge(UtcPeriod(recal_start, recal_stop)))
    set_calendar_formatter(Calendar())
    #plt.interactive(1)
    plt.title("Continuously recalibrated discharge vs observed")
    plt.xlabel("Time in UTC")
    plt.ylabel(r"Discharge in $\mathbf{m^3s^{-1}}$", verticalalignment="top", rotation="horizontal")
    plt.gca().yaxis.set_label_coords(0, 1.1)

 def plot_exon(row):
     start = int(row["exon_start"])
     stop = int(row["exon_stop"])
     size = stop - start    
     
     #print start, stop
     
     rectangle = plt.Rectangle((start, -20), size, 10, fc='red')
     plt.gca().add_patch(rectangle)

def build_normalized_histogram(digit_freq, labels):
    fig = plt.figure()
    plt.hist(digit_freq.keys(), weights=digit_freq.values())
    plt.title('Normalized Histogram')
    plt.xlabel('Digit #')
    plt.ylabel('Frequency')
    plt.gca().set_xlim([0, 9])
    fig.savefig('norm_hist.png')
    return

 def plot_column_elem_degree_distribution(self,logscale = True,bins = None):
     column_elem_degrees = self.degrees_of_column_elems()
     if bins is None:
         bins = int(numpy.max(column_elem_degrees))
     pylab.hist(column_elem_degrees,bins=bins)
     if logscale:
         pylab.gca().set_yscale("log")
         pylab.gca().set_xscale("log")
     pylab.show()