def draw_log_hist(X):
    X = X.tocsc().tocoo()  # collapse multiple records. I don't think it is needed

    # we are interested only in existence of a token in user posts, not it's quantity
    vf = np.vectorize(lambda x: 1 if x > 0 else 0)
    X_data_booleaned = vf(X.data)
    X = coo_matrix((X_data_booleaned, (X.row, X.col)), shape=X.shape)

    # now we will calculate (1, 1, ... 1) * X to sum up rows
    features_counts = np.ones(X.shape[0]) * X

    features_counts_sorted = np.sort(features_counts)
    features_counts_sorted = features_counts_sorted[::-1]  # this is how decreasing sort looks like in numpy
    ranks = np.arange(features_counts_sorted.size)

    plt.figure()
    plt.semilogy(ranks, features_counts_sorted,
                 color='red',
                 linewidth=2)
    plt.title('For each feature (word), how many users has it at least once?')
    plt.ylabel("number of users which has this word at least once")
    plt.xlabel("rank")
    # plt.show()

    return features_counts

def main():
    conn = krpc.connect()
    vessel = conn.space_center.active_vessel
    streams = init_streams(conn,vessel)
    print vessel.control.throttle
    plt.axis([0, 100, 0, .1])
    plt.ion()
    plt.show()

    t0 = time.time()
    timeSeries = []
    vessel.control.abort = False
    while not vessel.control.abort:

        t_now = time.time()-t0
        tel = Telemetry(streams,t_now)
        timeSeries.append(tel)
        timeSeriesRecent = timeSeries[-40:]

        plt.cla()
        plt.semilogy([tel.t for tel in timeSeriesRecent], [norm(tel.angular_velocity) for tel in timeSeriesRecent])
        #plt.semilogy([tel.t for tel in timeSeriesRecent[1:]], [quat_diff_test(t1,t2) for t1,t2 in zip(timeSeriesRecent,timeSeriesRecent[1:])])
        #plt.axis([t_now-6, t_now, 0, .1])
        plt.draw()
        plt.pause(0.0000001)
        #time.sleep(0.0001)

    with open('log.json','w') as f:
        f.write(json.dumps([tel.__dict__ for tel in timeSeries],indent=4))

    print 'The End'

def pal5():
    from streams import usys

    r0 = np.array([[8.312877511,0.242593717,16.811943627]])
    v0 = ([[-52.429087,-96.697363,-8.156130]]*u.km/u.s).decompose(usys).value

    x0 = np.append(r0,v0)
    acc = np.zeros((2,3))
    potential = LawMajewski2010()

    def F(t,X):
        x,y,z,px,py,pz = X.T
        dH_dq = potential._acceleration_at(X[...,:3], 2, acc)
        return np.hstack((np.array([px, py, pz]).T, dH_dq))

    integrator = DOPRI853Integrator(F)
    nsteps = 100000
    dt = 0.1
    print(nsteps*dt*u.Myr)
    nsteps_per_pullback = 10
    d0 = 1e-5

    LEs, xs = lyapunov(x0, integrator, dt, nsteps,
                       d0=d0, nsteps_per_pullback=nsteps_per_pullback)

    print("Lyapunov exponent computed")
    plt.clf()
    plt.semilogy(LEs, marker=None)
    plt.savefig('pal5_le.png')

    plt.clf()
    plt.plot(np.sqrt(xs[:,0]**2+xs[:,1]**2), xs[:,2], marker=None)
    plt.savefig('pal5_orbit.png')

def plot_magnif_from_time_terms(time_terms, magnifs, fmt="--ro", magnif_log_scale=False):
    plt.xlabel("time term, (t - t_max)/t_E")
    plt.ylabel("magnification")
    if magnif_log_scale:
        plt.semilogy(time_terms, magnifs, fmt)
    else:
        plt.plot(time_terms, magnifs, fmt)

	def __call__(self,u,v,w,iteration):
		q = 4

		plt.cool()
		if self.x == None:
			ny = v.shape[1]
			nz = v.shape[0]
			self.x,self.y = np.meshgrid(range(ny),range(nz))
		x,y = self.x,self.y

		if self.iterations == None:
			self.iterations = self.sim.bulk_calc(getIteration())
		all_itr = self.iterations

		if self.xvar == None:
			class temp(sim_operation):
				def get_params(self):
					return ["u"]
				def __call__(self,u):
					return np.max(self.sim.ddx(u))

			self.xvar = self.sim.bulk_calc(temp())
		xvar_series = self.xvar

		min = np.min(xvar_series)
		max = np.max(xvar_series)
		if min <= 0:
			min = 0.000001
		if max <= min:
			max = 0.00001
	
		avgu = np.average(u,2)
		avgv = np.average(v,2)
		avgw = -np.average(w,2)
		xd = self.sim.ddx(u)
		xd2d = np.max(xd,2)
		xd1d = np.max(xd2d,1)

		plt.subplot(221)
		plt.imshow(avgu)
		plt.quiver(x[::q,::q],y[::q,::q],avgv[::q,::q],avgw[::q,::q])
		plt.title('Avg u')
		plt.axis("tight")

		plt.subplot(222)
		plt.imshow(xd2d)
		plt.title('Max x Variation (y-z)')
		plt.axis("tight")

		plt.subplot(223)
		plt.plot(xd1d)
		plt.title('Max x Variation (z)')
		plt.axis("tight")

		plt.subplot(224)
		plt.plot(all_itr,xvar_series, '--')
		plt.plot([iteration,iteration],[min,max])
		plt.semilogy()
		plt.title('Max x Variation (t)')
		plt.axis("tight")

def plot_gens(images, rowlabels, losses):
    '''
    From great jupyter notebook by Tim Sainburg:
    http://github.com/timsainb/Tensorflow-MultiGPU-VAE-GAN
    '''
    examples = 8
    fig, ax = plt.subplots(nrows=len(images), ncols=examples, figsize=(18, 8))
    for i in range(examples):
        for j in range(len(images)):
            ax[(j, i)].imshow(create_image(images[j][i]), cmap=plt.cm.gray,
                              interpolation='nearest')
            ax[(j, i)].axis('off')
    title = ''
    for i in rowlabels:
        title += ' {}, '.format(i)
    fig.suptitle('Top to Bottom: {}'.format(title))
    plt.show()
    #fig.savefig(''.join(['imgs/test_',str(epoch).zfill(4),'.png']),dpi=100)
    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(20, 10), linewidth = 4)

    D_plt, = plt.semilogy((losses['discriminator']), linewidth=4, ls='-',
                          color='b', alpha=.5, label='D')
    G_plt, = plt.semilogy((losses['generator']), linewidth=4, ls='-',
                          color='k', alpha=.5, label='G')

    plt.gca()
    leg = plt.legend(handles=[D_plt, G_plt],
                     fontsize=20)
    leg.get_frame().set_alpha(0.5)
    plt.show()

def plot_trunc_gr_model(aval, bval, min_mag, max_mag, dmag, catalogue=None,
        completeness=None, figure_size=None, filename=None, filetype='png', 
        dpi=300):
    """
    Plots a Gutenberg-Richter model
    """
    input_model = TruncatedGRMFD(min_mag, max_mag, dmag, aval, bval)
    if not catalogue:
        # Plot only the modelled recurrence
        annual_rates, cumulative_rates = _get_recurrence_model(input_model)
        plt.semilogy(annual_rates[:, 0], annual_rates[:, 1], 'b-')
        plt.semilogy(annual_rates[:, 0], cumulative_rates, 'r-')
        plt.xlabel('Magnitude', fontsize='large')
        plt.ylabel('Annual Rate', fontsize='large')
        plt.legend(['Incremental Rate', 'Cumulative Rate'])
        _save_image(filename, filetype, dpi)
    else:
        completeness = _check_completeness_table(completeness, catalogue)
        plot_recurrence_model(input_model,
                              catalogue,
                              completeness,
                              input_model.bin_width,
                              figure_size,
                              filename,
                              filetype,
              		      dpi)

def plota_teste5(arqsaida):
    n, c, t = np.loadtxt(arqsaida, unpack=True)

    # Calcula os coeficientes de um ajuste a um polinômio de grau 2 usando
    # o método dos mínimos quadrados
    coefs = np.polyfit(n, c, 2)
    p = np.poly1d(coefs)

    # set_yscale('log')
    # set_yscale('log')
    plt.semilogy(n, p(n), label='$n^2$')
    plt.semilogy(n, c, 'ro', label='bubble sort')

    # Posiciona a legenda
    plt.legend(loc='upper left')

    # Posiciona o título
    plt.title('Análise da complexidade de \ntempo do método da bolha')

    # Rotula os eixos
    plt.xlabel('Tamanho do vetor (n)')
    plt.ylabel('Número de comparações')

    plt.savefig('bubble5.png')
    plt.show()

def plot_prob_sums(prob_sums, folder):
    """Plot the probability sums of the original motif against iterations."""
    # Create the folder if it does not exist
    if not os.path.exists("figures/%s" % folder):
        os.makedirs("figures/%s" % folder)

    # Plot on automatic axis
    figure_file = "figures/%s/probability_mass.png" % folder
    plt.title("%s | Prob Mass of Original Motif" % folder)
    plt.ylabel("Probability Mass")
    plt.xlabel("Iterations")
    plt.grid(True)
    plt.plot(prob_sums)
    pylab.savefig(figure_file, format="png")
    plt.clf()
    
    # Plot on semi-log axis
    figure_file = "figures/%s/probability_mass_semilog.png" % folder
    plt.title("%s | Prob Mass of Original Motif" % folder)
    plt.xlabel("Iterations")
    plt.ylabel("Probability Mass (semi-log)")
    plt.grid(True)
    plt.semilogy(prob_sums)
    pylab.savefig(figure_file, format="png")
    plt.clf()

def plot_KLs(KLs, folder, filename, title):
    if not os.path.exists("figures/%s" % folder):
        os.makedirs("figures/%s" % folder)

    # Plot on automatic axis
    figure_file = "figures/%s/%s.png" % (folder, filename)
    plt.title("%s | %s" % (folder, title))
    plt.ylabel("KL")
    plt.xlabel("Iterations")
    plt.grid(True)
    plt.plot(KLs)
    pylab.savefig(figure_file, format="png")
    plt.clf()
    
    # Plot on semi-log axis
    try:
        figure_file = "figures/%s/%s_semilog.png" % (folder, filename)
        plt.title("%s | %s" % (folder, title))
        plt.xlabel("Iterations")
        plt.ylabel("KL (semi-log)")
        plt.grid(True)
        plt.semilogy(KLs)
        pylab.savefig(figure_file, format="png")
    except:
        debug("Could not generate semilog plot.")
    plt.clf()

def on_off_experiment2(num_motifs=100,filename="gini-vs-mi-correlation-in-on-off-spoofs.pdf"):
    """compare MI vs Gini on biological_motifs"""
    bio_motifs = [getattr(tfdf,tf) for tf in tfdf.tfs]
    Ns = map(len, bio_motifs)
    spoofses = [spoof_on_off_motif(motif,num_motifs=num_motifs,trials=1) for motif in bio_motifs]
    spoof_ginises = mmap(motif_gini,tqdm(spoofses))
    spoof_mises = mmap(total_motif_mi,tqdm(spoofses))
    cors, ps = [],[]
    for ginis, mis in zip(ginises, mises):
        cor, p = pearsonr(ginis,mis)
        cors.append(cor)
        ps.append(p)
    q = fdr(ps)
    
    plt.scatter(cors,ps,filename="gini-vs-mi-correlation-in-on-off-spoofs.pdf")
    plt.plot([-1,1],[q,q],linestyle='--',label="FDR-Adjusted Significance Level")
    plt.semilogy()
    plt.legend()
    plt.xlabel("Pearson Correlation Coefficient")
    plt.ylabel("P value")
    plt.xlim([-1,1])
    plt.ylim([10**-4,1+1])
    cor_ps = zip(cors,ps)
    sig_negs = [(c,p) for (c,p) in cor_ps if c < 0 and p < q]
    sig_poses = [(c,p) for (c,p) in cor_ps if c > 0 and p < q]
    insigs = [(c,p) for (c,p) in cor_ps if p > q]
    def weighted_correlation(cor_p_Ns):
        cors,ps,Ns = transpose(cor_p_Ns)
        return sum([cor*N for (cor,N) in zip (cors,Ns)])/sum(Ns)
    plt.title("Gini-MI Correlation Coefficient vs. P-value for On-Off Simulations from Prokaryotic Motifs")
    maybesave(filename)

def plot_cmf_with_arbitrary_input(catalog, cmf_input, bins=20, hist_range=(4, 7), 
                                  mass_column_name='mass', **kwargs):

    fig = plt.figure()

    number_in_bin_q2, bin_edges, ch = plt.hist(np.log10(catalog[mass_column_name]), cumulative=-1, log=True, bins=bins, range=hist_range)
    bin_centers_q2 = (bin_edges[1:] + bin_edges[:-1])/2.
    plt.clf()
    plt.plot(bin_centers_q2, number_in_bin_q2, 'ko' )
    plt.semilogy()
    plt.xlabel(r"log$_{10}$ (M$_{GMC}$ / M$_\odot$)")
    plt.ylabel("n(M > M')")

    M_0, N_0, gamma = cmf_input

    m_array = np.linspace(min(bin_edges), max(bin_edges), 50)
    n_array = truncated_cloudmass_function([M_0, N_0, gamma], 10**m_array)

    plt.plot(m_array, n_array, label="$\\gamma = {0:.2f}$,\n$M_0={1:.2e}$,\n$N_0={2:.1f}$".format(gamma, M_0, N_0))

    text_string = r"$N(M' > M) = N_0 \left [ \left ( \frac{M}{M_0} \right )^{\gamma+1} - 1 \right ]$"

    plt.text(4.1, 3, text_string, fontsize=18)
    plt.xlim(*hist_range)
    plt.ylim(0.7, 1e3)

    plt.legend(loc='upper right')

    return fig    

def throughputs(output='throughputs.pdf'):
    """
    Plot throughputs, compares to HST WFC3 UVIS F600LP

    :param output: name of the output file
    :type output: str

    :return: None
    """
    #comparison
    bp1 = S.ObsBandpass('wfc3,uvis2,f600lp')

    #VIS
    bp, bpEoL = _VISbandpass()

    #ghost
    bpG, bpEoLG = _VISbandpassGhost()

    #plot
    plt.semilogy(bp1.wave/10., bp1.throughput, 'r-', label='WFC3 F600LP')
    plt.semilogy(bp.wave/10., bp.throughput, 'b-', label='VIS Best Estimate')
    plt.semilogy(bpEoL.wave/10., bpEoL.throughput, 'g--', label='VIS EoL Req.')
    plt.semilogy(bpG.wave/10., bpG.throughput, 'm-', label='VIS Ghost')
    plt.semilogy(bpEoLG.wave/10., bpEoLG.throughput, 'y-.', label='VIS Ghost EoL')
    plt.xlim(230, 1100)
    plt.xlabel(r'Wavelength [nm]')
    plt.ylabel(r'Total System Throughput')
    plt.legend(shadow=True, fancybox=True, loc='best')
    plt.savefig(output)
    plt.close()

def test_airy_1d(display=False):
    """ Compare analytic airy function results to the expected locations
    for the first three dark rings and the FWHM of the PSF."""
    lam = 1.0e-6
    D = 1.0
    r, airyprofile = airy_1d(wavelength=lam, diameter=D, length=20480, pixelscale=0.0001)

    # convert to units of lambda/D
    r_norm = r*_ARCSECtoRAD / (lam/D)
    if display:
        plt.semilogy(r_norm,airyprofile)
        plt.axvline(1.028/2, color='k', ls=':')
        plt.axhline(0.5, color='k', ls=':')
        plt.ylabel('Intensity relative to peak')
        plt.xlabel('Separation in $\lambda/D$')
        for rad in airy_zeros:
            plt.axvline(rad, color='red', ls='--')

    airyfn = scipy.interpolate.interp1d(r_norm, airyprofile)
    # test FWHM occurs at 1.028 lam/D, i.e. HWHM is at 0.514
    assert (airyfn(0.5144938) - 0.5) < 1e-5

    # test first minima occur near 1.22 lam/D, 2.23, 3.24 lam/D
    # TODO investigate/improve numerical precision here?
    for rad in airy_zeros:
        #print(rad, airyfn(rad), airyfn(rad+0.005))
        assert airyfn(rad) < airyfn(rad+0.0003)
        assert airyfn(rad) < airyfn(rad-0.0003)

def create_plot(freq_mat, data_mat, length_vec, color_keys, plot_name, output_dir=""):
	print("Creating plot: {0:s}".format(plot_name) )
	pl.figure(1)
	pl.clf()
	pl.hold(True)
	
	for idx, freq in enumerate(freq_mat):
		data = data_mat[idx]
		length = length_vec[idx]
		
		num_pos = 0
		num_neg = 0
		for el in data:
			if el > 0:
				num_pos += 1
			if el < 0:
				num_neg +=1
		if num_pos > 0:
			pl.semilogy(freq/1e9, data, color_keys[length])
		if num_neg > 0:
			pl.semilogy(freq/1e9, -data, color_keys[length] + ":")
		pl.xlabel('Frequency (GHz)')
		pl.ylabel("PUL R ($\Omega$/m)")
		
	output_name = os.path.join(output_dir, plot_name)
	pl.savefig(output_name)

def do_halomass_plots(fosc):
    """Plot halo mass, distance to the halo and the relationship between eq. width and column density"""
    ahalos = {
                    4:CIVPlottingSpectra(5, myname.get_name(4, box=25), None, None, savefile="rand_civ_spectra.hdf5", spec_res=5.,label=labels[4]+" 25"),
                   9:CIVPlottingSpectra(5, myname.get_name(9, box=25), None, None, savefile="rand_civ_spectra.hdf5", spec_res=5.,label=labels[9]+" 25"),
                    7:CIVPlottingSpectra(5, myname.get_name(7, box=25), None, None, savefile="rand_civ_spectra.hdf5", spec_res=5.,label=labels[7]+" 25")}
    for (ll, ahalo) in ahalos.items():
        ahalo.plot_eqw_mass("C",4,1548,color=colors[ll])
    plt.legend(loc="upper left")
    save_figure(path.join(outdir,"civ_eqwvsmass"))
    plt.clf()
    for (ll, ahalo) in ahalos.items():
        ahalo.plot_eqw_dist("C",4,1548,color=colors[ll])
    plt.legend(loc="upper left")
    save_figure(path.join(outdir,"civ_eqwvsdist"))
    plt.clf()
    ahalos['I']=CIVPlottingSpectra(68, path.expanduser("~/data/Illustris"), None, None, savefile="rand_civ_spectra.hdf5", spec_res=5.,label=labels['I']+" 75")
    ccc = {1e12: "yellow", 1e15:"red"}
    for tag in ('I', 4, 7):
        for nmin in (1e12, 1e15):
            nminstr = str(np.round(np.log10(nmin),1))
            #for (ll, ahalo) in ahalos.iteritems():
            ahalos[tag].plot_mass_hist(elem="C",ion=4,color=ccc[nmin], ls=lss[tag],nmin=nmin, label=ahalos[tag].label+" "+nminstr)
    plt.legend(loc="upper left",ncol=1)
    save_figure(path.join(outdir,"civ_halos_hist"))
    plt.clf()
    ahalos['I'].plot_eq_width_vs_col_den("C",4,1548)
    eqw = np.linspace(-3, 0.5,50)
    plt.semilogy(eqw, linear_cog_col(10**eqw, 1548, fosc), '-',color="black")
    plt.ylim(1e12,1e16)
    plt.xlim(-2.5,0.5)
    save_figure(path.join(outdir,"civ_eqwvscolden"))
    plt.clf()

def runAlgorithm(data, d, k):
    W = None

    alphas = 10.0 ** linspace(-4, 1, 6)
    #alphas = 10.0 ** array([-2,-1,0,-2,-3,-2])
    expNum = len(alphas)
    allStates = []
    for (alpha, i) in zip(alphas, count(1)):
        states = jointlyPenalizedMultiplePCA(data, d, alpha, W, k=k)
        allStates.append(states)
        weightVectors = [s.weights for s in states]
        W = array(weightVectors)

    figure()
    for (states, alpha, i) in zip(allStates, alphas, count(1)):
        subplot(expNum, 1, i)
        weightVectors = [s.weights for s in states]
        W = array(weightVectors)
        plot(W.T)
        title(r'Run with $\alpha $ = %f' % alpha)

    figure()
    for (states, alpha, i) in zip(allStates, alphas, count(1)):
        subplot(expNum, 1, i)
        lossVectors = [s.squaredLosses() for s in states]
        L = array(lossVectors)
        semilogy(L.T)
        title(r'Run with $\alpha $ = %f' % alpha)

    show()

 def advance(self, t, plotresult=False):
     y0 = self.concs * self.molWeight
     y0 = append(y0, self.thickness)
     yt = odeint(self.rightSideofODE, y0, t)
     if (plotresult):
         import matplotlib.pyplot as plt
         plt.figure()
         plt.axes([0.1, 0.1, 0.6, 0.85])
         plt.semilogy(t, yt)
         plt.ylabel('mass concentrations (kg/m3)')
         plt.xlabel('time(s)')
         #plt.legend(self.speciesnames)
         for i in range(len(self.speciesnames)):
             plt.annotate(
                 self.speciesnames[i], (t[-1], yt[-1, i]),
                 xytext=(20, -5),
                 textcoords='offset points',
                 arrowprops=dict(arrowstyle="-"))
         plt.show()
     self.thickness = yt[-1][-1]
     ytt = yt[-1][:-1]
     #        for iii in range(len(ytt)):
     #            if ytt[iii]<0:
     #                ytt[iii]=0.
     molDens = ytt / self.molWeight
     self.concs = molDens
     self.molFrac = molDens / sum(molDens)
     self.massFrac = ytt / sum(ytt)

 def plot(self):
     from matplotlib import pyplot
     pyplot.figure()
     it_ls = [] #ls = lineseach
     y_ls = []
     it_ga = [] #gradient approximation
     y_ga = []
     gradApprox = False
     for i in range(len(self.x)):
         y = norm( self.f_x[i] ) + 10**-9
         if self.notes.has_key(i):
             if self.notes[i] == 'starting gradient approximation':
                 gradApprox = True
             if self.notes[i] == 'finished gradient approximation':
                 gradApprox = False
         if gradApprox:
             it_ga.append( i )
             y_ga.append( y ) 
         else:
             it_ls.append( i )
             y_ls.append( y )
     pyplot.semilogy( it_ls, y_ls, 'go') 
     pyplot.semilogy( it_ga, y_ga, 'bx') 
     pyplot.xlabel('function evaluation')
     pyplot.ylabel('norm(f(x)) + 10**-9')
     pyplot.legend(['line searches', 'gradient approx' ])
                   
     pyplot.show()

def compute_lyapunov(orbit_name, halo_shape='oblate'):
    if halo_shape == "prolate":
        params = prolate_params
        raise NotImplementedError()
    elif halo_shape == "oblate":
        params = oblate_params
    else:
        raise ValueError("Invalid halo shape.")

    integrator = DOPRI853Integrator(F, func_args=params)
    nsteps = 10000
    dt = 1.
    print(nsteps*dt*u.Myr)
    nsteps_per_pullback = 10
    d0 = 1e-5

    LEs, xs = lyapunov(x0, integrator, dt, nsteps,
                       d0=d0, nsteps_per_pullback=nsteps_per_pullback)

    print("Lyapunov exponent computed")
    plt.figure(figsize=(10,10))
    plt.clf()
    plt.semilogy(LEs, marker=None)
    plt.savefig('zotos_le_{}.png'.format(orbit_type))

    plt.clf()
    plt.plot(xs[:,0], xs[:,1], marker=None)
    plt.xlim(0,15)
    plt.ylim(-15,15)
    plt.savefig('zotos_orbit_{}.png'.format(orbit_type))