def show_ratios(cpu):
    cpu.ratios.sort(key=lambda x: x[0])
    pl.figure("Tuning")
    pl.plot([x[0] for x in cpu.ratios], [x[2] for x in cpu.ratios])
    
    pl.figure("Tuning samples")
    pl.scatter([x[0] for x in cpu.ratios], [x[2] * x[0] for x in cpu.ratios])

def fit(w, f, e, mw, mf, vgrid, npol,
        sigrange=None, vrange=None, doppler=doppler, plot=False):

    vgrid = Quantity(vgrid, u.km/u.s)
    chi2 = Table([vgrid.value, np.zeros_like(vgrid.value)], names=['v','chi2'])
    chi2['v'].units = vgrid.unit

    fit1 = Fit1(w, f, e, mw, mf, npol, doppler)

    chi2['chi2'] = np.array([fit1(v)[0] for v in vgrid])

    chi2.meta['ndata'] = len(f)
    chi2.meta['npar'] = npol+1+1
    chi2.meta['ndof'] = chi2.meta['ndata']-chi2.meta['npar']

    if plot:
        import matplotlib.pylab as plt
        plt.scatter(chi2['v'], chi2['chi2'])

    if vrange is None and sigrange is None or len(vgrid) < 3:
        ibest = chi2['chi2'].argmin()
        vbest, bestchi2 = chi2[ibest]
        chi2.meta['vbest'] = vbest
        chi2.meta['verr'] = 0.
        chi2.meta['bestchi2'] = bestchi2
    else:
        vbest, verr, bestchi2 = minchi2(chi2, vrange, sigrange, plot=plot)

    _, fit, mfi = fit1(vbest)
    chi2.meta['wmean'] = fit1.wmean
    chi2.meta['continuum'] = fit1.sol
    return chi2, fit, mfi

    def visualization2(self, sp_to_vis=None):
        if sp_to_vis:
            species_ready = list(set(sp_to_vis).intersection(self.all_sp_signatures.keys()))
        else:
            raise Exception('list of driver species must be defined')

        if not species_ready:
            raise Exception('None of the input species is a driver')

        for sp in species_ready:
            # Setting up figure
            plt.figure()
            plt.subplot(313)

            mon_val = OrderedDict()
            signature = self.all_sp_signatures[sp]
            for idx, mon in enumerate(list(set(signature))):
                if mon[0] == 'C':
                    mon_val[self.all_comb[sp][mon] + (-1,)] = idx
                else:
                    mon_val[self.all_comb[sp][mon]] = idx

            mon_rep = [0] * len(signature)
            for i, m in enumerate(signature):
                if m[0] == 'C':
                    mon_rep[i] = mon_val[self.all_comb[sp][m] + (-1,)]
                else:
                    mon_rep[i] = mon_val[self.all_comb[sp][m]]
            # mon_rep = [mon_val[self.all_comb[sp][m]] for m in signature]

            y_pos = numpy.arange(len(mon_val.keys()))
            plt.scatter(self.tspan[1:], mon_rep)
            plt.yticks(y_pos, mon_val.keys())
            plt.ylabel('Monomials', fontsize=16)
            plt.xlabel('Time(s)', fontsize=16)
            plt.xlim(0, self.tspan[-1])
            plt.ylim(0, max(y_pos))

            plt.subplot(312)

            for name in self.model.odes[sp].as_coefficients_dict():
                mon = name
                mon = mon.subs(self.param_values)
                var_to_study = [atom for atom in mon.atoms(sympy.Symbol)]
                arg_f1 = [numpy.maximum(self.mach_eps, self.y[str(va)][1:]) for va in var_to_study]
                f1 = sympy.lambdify(var_to_study, mon)
                mon_values = f1(*arg_f1)
                mon_name = str(name).partition('__')[2]
                plt.plot(self.tspan[1:], mon_values, label=mon_name)
            plt.ylabel('Rate(m/sec)', fontsize=16)
            plt.legend(bbox_to_anchor=(-0.1, 0.85), loc='upper right', ncol=1)

            plt.subplot(311)
            plt.plot(self.tspan[1:], self.y['__s%d' % sp][1:], label=parse_name(self.model.species[sp]))
            plt.ylabel('Molecules', fontsize=16)
            plt.legend(bbox_to_anchor=(-0.15, 0.85), loc='upper right', ncol=1)
            plt.suptitle('Tropicalization' + ' ' + str(self.model.species[sp]))

            # plt.show()
            plt.savefig('s%d' % sp + '.png', bbox_inches='tight', dpi=400)

def scipy_stuff():
  from scipy.interpolate import griddata
  from matplotlib import pylab
  import cPickle as pickle
  print "loading points"
  points, x_diff, y_diff = pickle.load(open("temp_data.pickle", "rb"))

  y_pts, x_pts = zip(*points)

  print "Creating grid points"
  grid_points = []
  for j in range(2500):
    for i in range(2500):
      grid_points.append((j, i))

  print "Gridding data"
  x_grid = griddata(points, x_diff, grid_points)
  y_grid = griddata(points, y_diff, grid_points)
  x_grid.shape = (2500, 2500)
  y_grid.shape = (2500, 2500)

  print "Plotting"
  pylab.subplot(3, 1, 1)
  pylab.imshow(x_grid)
  pylab.subplot(3, 1, 2)
  pylab.imshow(y_grid)
  pylab.subplot(3, 1, 3)
  pylab.scatter(x_pts, y_pts)
  pylab.show()

def fit_plot_unlabeled_data(unlabeled_data_x, labeled_data_x, labeled_data_y, fit_order, data_type, other_data_list, other_data_name):
    output = open('predictions.csv','wb')
    coeffs = np.polyfit(labeled_data_x, labeled_data_y, fit_order) #does poly git to nth deg on labeled data
    fit_eq = np.poly1d(coeffs) #Eqn from fit
    predicted_y = fit_eq(unlabeled_data_x)
    i = 0
    writer = csv.writer(output,delimiter=',')
    header = [str(data_type),str(other_data_name),'Predicted_Num_Inc']
    writer.writerow(header)
    while i < len(predicted_y):
        output_data = [unlabeled_data_x[i],other_data_list[i],predicted_y[i]]
        writer.writerow(output_data)
        print 'For '+str(data_type)+' of: '+str(unlabeled_data_x[i])+', Predicted Number of Incidents is: '+str(predicted_y[i])
        i = i + 1
    plt.scatter(unlabeled_data_x, predicted_y, color='blue', label='Predicted Number of Incidents')
    fit_line_x = np.arange(min(unlabeled_data_x), max(unlabeled_data_x), 1)
    plt.plot(fit_line_x, fit_eq(fit_line_x), color='red',linestyle='dashed',label=' Order '+str(fit_order)+' Polynomial Fit')
#____Use below line to plot actual data also!! 
    #plt.scatter(labeled_data_x, labeled_data_y, color='green', label='Actual Incident Report Data')
    plt.title('Predicted Number of 311 Incidents by '+str(data_type))
    plt.xlabel(str(data_type))
    plt.ylabel('Number of 311 Incidents')
    plt.grid()
    plt.xlim([min(unlabeled_data_x)-1500, max(unlabeled_data_x)+1500])
    plt.legend(loc='upper left')
    plt.show()

def plot_prediction_accuracy(x, y):
    plt.scatter(x, y, c='g', alpha=0.5)
    plt.title('Logistic Regression')
    plt.xlabel('r')
    plt.ylabel('Prediction Accuracy')
    plt.xlim(0,200)
    plt.show()

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 main():
    # an example non-autonomous function
    x0 = 1
    t = np.linspace(1,3,500)

    # use the same syntax as odeint
    sol = LSolve(example,x0,t,args=(1,1))
    
    if matplotlib_module:
        mp.figure(1)
        mp.title("Example solution")
        mp.plot(t,sol)


    # example integrate and fire code
    x0 = 0
    t2 = np.linspace(0,10,500)
    
    # again the syntax is the same as odeint, but we add aditional inputs,
    # including a flag to track spikes (IF models only):
    threshold = 5
    sol2,spikes = LSolve(QIF,x0,t,threshold=threshold,reset=0,spike_tracking=True,args=(5,))

    # extract spike times
    spikes[spikes==0.0] = None
    spikes[spikes==1.0] = threshold

    if matplotlib_module:
        mp.figure(2)
        mp.title("QIF model with noise")
        mp.plot(t2,sol2)
        mp.scatter(t2,spikes,color='red',facecolor='red')
        mp.show()

def movie_plotter(components, movies, movie_id="all", x_buffer=3, y_buffer=2):
    if movie_id == "all":
        plt.scatter(components[:,0], components[:,1])
        plt.xlabel("Component 1")
        plt.ylabel("Component 2")
        plt.show()
    else:
        x = components[movie_id][0]
        y = components[movie_id][1]

        xs = [x - x_buffer, x + x_buffer]
        ys = [y - y_buffer, y + y_buffer]

        plt.scatter(components[:,0], components[:,1])
        plt.xlim(xs)
        plt.ylim(ys)
        plt.xlabel("Component 1")
        plt.ylabel("Component 2")

        for x, y, title in zip(components[:,0], components[:,1], movies['movie_title']):
            if x >= xs[0] and x <= xs[1] and y >= ys[0] and y <= ys[1]:
                try:
                    plt.text(x, y, title)
                except:
                    pass

def classification_regions(network, title, img_file_name, interval=100):
    coords = [
        (i / interval, j / interval)
        for j
        in range(0, interval + 1, 1)
        for i
        in range(0, interval + 1, 1)]

    classified_records = []

    for coord in coords:
        output = network.run(coord)

        classified_records.append(
            [coord, output.index(max(output)) + 1])

    plt.scatter(
        [record[0][0] for record in classified_records],
        [record[0][1] for record in classified_records],
        c=[record[1]    for record in classified_records],
    )

    plt.xlim((0, 1))
    plt.ylim((0, 1))

    plt.title(title)
    plt.xlabel('Six-fold rotational symmetry')
    plt.ylabel('Eccentricity')

    plt.savefig(img_file_name)

def main_k_nearest_neighbour(k):
    X, y = make_blobs(n_samples=100,
                      n_features=2,
                      centers=2,
                      cluster_std=1.0,
                      center_box=(-10.0, 10.0))

    h = .4
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))

    z = np.c_[xx.ravel(), yy.ravel()]

    z_f = []
    for i_z in z:
        z_f.append(k_nearest_neighbour(X, y, i_z, k, False))

    zz = np.array(z_f).reshape(xx.shape)

    plt.figure()
    plt.contourf(xx, yy, zz, cmap=plt.cm.Paired)
    plt.axis('tight')
    plt.scatter(X[:, 0], X[:, 1], c=y)

    plt.show()

def winding_number(filename,window,bins=0):
    
    data=tools.read_matrix_from_file(filename);
    print "file read."
    print len(data[1])

    if (bins ==0):
        bins=window
    
    times=np.zeros(bins)
    values=np.zeros(bins)
    ns=np.zeros(bins)
    step=window/bins
    
    for i in range(0,bins):
        times[i]=i*step
    
    for k in range(0,len(data[1])-window,window):
    
        for j in range(k,window+k,step):
        
            for i in range(0,bins):
                
                values[i]=values[i]+(data[1][j]-data[1][j - i*step])**2
                ns[i]=ns[i]+1
    
    for i in range(0,bins):
        if (ns[i] != 0):
            values[i]=values[i]/ns[i]
        else:
            values[i]=0
    plt.scatter(times,values)
    return [times,values]

 def plot_pulses(results, ymin=0, ymax=20):
     plt.plot(results["times"], results["amounts"])
     s = np.array([1] * len(results["times"]))
     c = np.array(["k"] * len(results["times"]))
     if "durations" in results:
         # semi-Markovian
         start = 0
         for d, pulse in zip(results["durations"],
                             results["pulses"]):
             end = min(start + d, len(results["times"]) - 1)
             if pulse:
                 c[start:end] = "red"
                 s[start:end] = 2
             start += d
     else:
         # Markovian
         for n, t in enumerate(results["times"]):
             pulse = results["pulses"][n]
             if pulse:
                 c[n] = "red"
                 s[n] = 2
     plt.scatter(results["times"], [1] * len(results["times"]), color=c, s=s)
     plt.xlabel(r"Time, $t$")
     plt.ylabel("Glucose amount")
     plt.ylim([ymin, ymax])
     plt.xlim([time_obj.t.min(), time_obj.t.max()])
     sns.despine()

def main():   
    x0 = np.loadtxt('ex/ex5Linx.dat')
    y = np.loadtxt('ex/ex5Liny.dat')
    x0.shape=x0.size,1
    y.shape = y.size,1
    
    plt.scatter(x0,y)
    x = polynomial_linear(x0)
#    x,mns,sstd = z_scale(x)

    theta_normal = linear_normal_equation(x,y, 1.0)
    print 'normal equation:'
    print theta_normal
    plot_fitting(theta_normal)
    plt.show()  
    m,n=x.shape 
    alphas = ( 0.01, 0.03, 0.1, 0.3, 1, 1.3 )    # if alpha >=1.3, no convergence result
    lambdas = (0, 1, 10)
    MAX_ITR = 100 
    for lam in lambdas:
        for alpha in alphas:
            theta,Js = linear_regression(x,y, MAX_ITR, alpha, lam)    
            if alpha==0.03 and lam==1:
                theta_best = theta
            plt.plot(Js)
        plt.xlabel('iterations')
        plt.ylabel('cost: J')
        plt.legend(['alpha: %s' %i for i in alphas])
        plt.show()

    print 'best theta in alpha:\n ', theta_best
    test = x0[-1]
    test.shape=test.size,1
    test = polynomial_linear(test)
    print 'predict of %s is %s' %(test, predict_linear(theta, test))

 def scatter(title, file_name, x_array, y_array, size_array, x_label, \
             y_label, x_range, y_range, print_pdf):
     '''
     Plots the given x value array and y value array with the specified 
     title and saves with the specified file name. The size of points on
     the map are proportional to the values given in size_array. If 
     print_pdf value is 1, the image is also written to pdf file. 
     Otherwise it is only written to png file.
     '''
     rc('text', usetex=True)
     rc('font', family='serif')
     plt.clf() # clear the ploting window, a must.                               
     plt.scatter(x_array, y_array, s =  size_array, c = 'b', marker = 'o', alpha = 0.4)
     if x_label != None:   
         plt.xlabel(x_label)
     if y_label != None:
         plt.ylabel(y_label)                
     plt.axis ([0, x_range, 0, y_range])
     plt.grid(True)
     plt.suptitle(title)
 
     Plotter.print_to_png(plt, file_name)
     
     if print_pdf:
         Plotter.print_to_pdf(plt, file_name)

def plot_trajectory(mu_vector):
    data0 = mu_vector[:, 0]
    data1 = mu_vector[:, 1]
    labels = ["{0}".format(i) for i in xrange(len(mu_vector))]
    plt.scatter(data0[:, 0], data0[:, 1], color="red")
    plt.scatter(data1[:, 0], data1[:, 1], color="blue")
    for i in xrange(len(mu_vector)):
        plt.annotate(
            labels[i],
            (data0[i, 0], data0[i, 1]),
            fontsize=5,
            xytext=(-10, 20),
            textcoords="offset points",
            ha="right",
            va="bottom",
            arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=0"),
        )
        plt.annotate(
            labels[i],
            (data1[i, 0], data1[i, 1]),
            fontsize=5,
            xytext=(-10, 20),
            textcoords="offset points",
            ha="right",
            va="bottom",
            arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=0"),
        )
    plt.savefig("Mean_Trajectory.png")
    plt.show()

def Decision_Surface(data, target, model, surface=True, probabilities=False, cell_size=.01):
    # Get bounds
    x_min, x_max = data[data.columns[0]].min(), data[data.columns[0]].max()
    y_min, y_max = data[data.columns[1]].min(), data[data.columns[1]].max()
    
    # Create a mesh
    xx, yy = np.meshgrid(np.arange(x_min, x_max, cell_size), np.arange(y_min, y_max, cell_size))
    meshed_data = pd.DataFrame(np.c_[xx.ravel(), yy.ravel()])
    
    # Add interactions
    for i in range(data.shape[1]):
        if i <= 1:
            continue
        
        meshed_data = np.c_[meshed_data, np.power(xx.ravel(), i)]

    if model != None:
        # Predict on the mesh
        if probabilities:
            Z = model.predict_proba(meshed_data)[:, 1].reshape(xx.shape)
        else:
            Z = model.predict(meshed_data).reshape(xx.shape)
    
    # Plot mesh and data
    if data.shape[1] > 2:
        plt.title("humor^(" + str(range(1,data.shape[1])) + ") and number_pets")
    else:
        plt.title("humor and number_pets")
    plt.xlabel("humor")
    plt.ylabel("number_pets")
    if surface and model != None:
        cs = plt.contourf(xx, yy, Z, cmap=plt.cm.coolwarm, alpha=0.4)
    color = ["blue" if t == 0 else "red" for t in target]
    plt.scatter(data[data.columns[0]], data[data.columns[1]], color=color)

def test(args):
    data = multivariate_normal([0, 0], [[1, 2], [2, 5]], int(args[1]))
    print(data)
    # PCA
    result = pca(data, base_num=int(args[2]))
    pc_base = result[0]
    print(pc_base)

    # Plotting
    fig = plt.figure()
    fig.add_subplot(1, 1, 1)
    plt.axvline(x=0, color="#000000")
    plt.axhline(y=0, color="#000000")
    # Plot data
    plt.scatter(data[:, 0], data[:, 1])
    # Draw the 1st principal axis
    pc_line = sp.array([-3.0, 3.0]) * (pc_base[1] / pc_base[0])
    plt.arrow(0, 0, -pc_base[0] * 2, -pc_base[1] * 2, fc="r", width=0.15, head_width=0.45)
    plt.plot([-3, 3], pc_line, "r")
    # Settings
    plt.xticks(size=15)
    plt.yticks(size=15)
    plt.xlim([-3, 3])
    plt.tight_layout()
    plt.show()
    plt.savefig("image.png")

    return 0

def plot_approx_error(r, error):
    plt.scatter(r, error, c='g', alpha=0.5)
    plt.title('Rank r Approximation Error')
    plt.xlabel('r')
    plt.ylabel('Frobenius Norm')
    plt.show()
    return

def plot_heatmap(x, y_mat, alpha=1,title=None,sizes=None,share_axis=False):
    if sizes is None:
        sizes = 60
    if y_mat.ndim == 1:
        y_mat = np.expand_dims(y_mat, 1)
    pl.close()
    fig = pl.figure(4)
    fig.suptitle(title)
    for index, y in enumerate(y_mat.T):
        if index == 0:
            ax1 = pl.subplot(y_mat.shape[1], 1, index + 1)
        else:
            if share_axis:
                pl.subplot(y_mat.shape[1], 1, index + 1, sharex=ax1, sharey=ax1)
            else:
                pl.subplot(y_mat.shape[1], 1, index + 1)
        red_values = normalize(y)
        I = np.isfinite(y) & np.isfinite(x[:,0]) & np.isfinite(x[:,1])
        colors = np.zeros((red_values.size, 3))
        colors[:,0] = red_values
        pl.ylabel(str(index))
        if I.mean > 0:
            print 'Percent skipped due to nans: ' + str(1-I.mean())
        pl.scatter(x[I,0], x[I,1], alpha=alpha, c=colors[I,:], edgecolors='none', s=sizes)
    move_fig(fig, 1000, 1000)
    pl.show(block=True)
    pass