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 get_OH(OIII4363,OIII4959,OIII5007,Hb):
	Te = np.arange(5000,20000,1)
	t3 = Te/1e4
	ne = 100	# cm^-3
	x = 1e-4*ne*t3**(-0.5)
	C_T = (8.44-1.09*t3+0.5*t3**2.-0.08*t3**3.)*(1.+0.0004*x)/(1.+0.044*x)
	log_OIII_ratio = 1.432/t3+np.log10(C_T)
	log_OIII_ratio_obs = np.log10((OIII4959+OIII5007)/OIII4363)

	Te_obs = []
	plt.clf()
	plt.plot(Te,log_OIII_ratio,color='black',marker='.',linestyle='none')
	plt.yscale('log')
	for i in range(len(OIII4363)):
		plt.axhline(log_OIII_ratio_obs[i],linestyle='--')
		d_ratio = abs(log_OIII_ratio_obs[i]-log_OIII_ratio)
		min_d_ratio = min(d_ratio)
		min_sub = list(d_ratio).index(min_d_ratio)
		Te_obs.append(Te[min_sub])
	plt.xlim(10000,20000)
	plt.ylim(1.,3)
	plt.xlabel('Te')
	plt.ylabel('(OIII4959+5007)/OIII4363')

	Te_obs = np.array(Te_obs)
	t3_obs = Te_obs/1e4
	logOIIIH = np.log10((OIII4959+OIII5007)/Hb)+6.200+1.251+1.251/t3_obs - \
				5*np.log10(t3_obs)-0.014*t3_obs
	t2_obs = -0.577+t3*(2.065-0.498*t3)
	logOIIH = np.log10(OII3727/Hb)+5.961+1.676/t2_obs-0.4*np.log10(t2_obs) - \
				0.034*t2_obs+np.log10(1+1.35*x)
	OH = 10**(logOIIIH-12.)+10**(logOIIIH-12.)
	logOH = 12 + np.log10(OH)

	return Te_obs,logOIIH,logOIIIH,logOH

def plot_ra(s1, s2, idxs=None, epsilon=0.25, fig=None):
    """Computes the RA plot of two groups of samples"""

    ## compute log2 values
    l1 = np.log2(s1 + epsilon)
    l2 = np.log2(s2 + epsilon)

    ## compute A and R
    r = l1 - l2
    a = (l1 + l2) * 0.5

    fig = pl.figure() if fig is None else fig
    pl.figure(fig.number)

    if idxs is None:
        pl.plot(a, r, '.k', markersize=2)
    else:
        pl.plot(a[~idxs], r[~idxs], '.k', markersize=2)
        pl.plot(a[idxs], r[idxs], '.r')

    pl.axhline(0, linestyle='--', color='k')

    pl.xlabel('(log2 sample1 + log2 sample2) / 2')
    pl.ylabel('log2 sample1 - log2 sample2')

    pl.tight_layout()

def testData(symbols, periods, low, high, portfolio, tradeSize, tradeFee, firstOnly=0, reverse=True, buyOnly=True, plot=False):     
    numSym = len(symbols)
    ress = [[] for sym in symbols]
    for i, sym in enumerate(symbols):
        data = pandas.read_csv(sym+".csv")
        if reverse:
            data = data[::-1]
        if firstOnly>0:
            res = data[:firstOnly].copy()
        else:
            res = data.copy()
        ress[i] = res
        rsiCompute(res, periods)
        buyOrSell(res, low, high)
    tp, nt, nb, ns = computeProfit(ress, symbols, portfolio, tradeSize, tradeFee, buyOnly)
    print("tp = ",tp,", nt = ",nt,", nb = ",nb,", ns = ",ns)
    print("profit % = ", tp/portfolio*100)

#     print(res['rsi'])
    if plot and numSym == 1:
        res = ress[0]
        n = len(res)
        rcParams['figure.figsize'] = 60, 12
        bs = np.array(res['bs'])
        x = np.arange(n)
        markers_onB = x[bs=='B']
        markers_onS = x[bs=='S']
        plt.plot(x, res["Adj Close"], '-bo', label='spy', markevery=markers_onB)
        plt.plot(x, res["rsi"], '-ro', label='RSI', markevery=markers_onS)
        plt.axhline(low)
        plt.axhline(high)
        plt.xlabel('Date')
        plt.ylabel('Price')
        plt.show()

def fancy_dendrogram(*args, **kwargs):
    '''
    Source: https://joernhees.de/blog/2015/08/26/scipy-hierarchical-clustering-and-dendrogram-tutorial/
    '''
    from scipy.cluster import hierarchy
    import matplotlib.pylab as plt
    
    max_d = kwargs.pop('max_d', None)
    if max_d and 'color_threshold' not in kwargs:
        kwargs['color_threshold'] = max_d
    annotate_above = kwargs.pop('annotate_above', 0)

    ddata = hierarchy.dendrogram(*args, **kwargs)

    if not kwargs.get('no_plot', False):
        plt.title('Hierarchical Clustering Dendrogram (truncated)')
        plt.xlabel('sample index or (cluster size)')
        plt.ylabel('distance')
        for i, d, c in zip(ddata['icoord'], ddata['dcoord'], ddata['color_list']):
            x = 0.5 * sum(i[1:3])
            y = d[1]
            if y > annotate_above:
                plt.plot(x, y, 'o', c=c)
                plt.annotate("%.3g" % y, (x, y), xytext=(0, -5),
                             textcoords='offset points',
                             va='top', ha='center')
        if max_d:
            plt.axhline(y=max_d, c='k')
    return ddata

def showKernel(dataOrMatrix, fileName = None, useLabels = True, **args) :
 
    labels = None
    if hasattr(dataOrMatrix, 'type') and dataOrMatrix.type == 'dataset' :
	data = dataOrMatrix
	k = data.getKernelMatrix()
	labels = data.labels
    else :
	k = dataOrMatrix
	if 'labels' in args :
	    labels = args['labels']

    import matplotlib

    if fileName is not None and fileName.find('.eps') > 0 :
        matplotlib.use('PS')
    from matplotlib import pylab

    pylab.matshow(k)
    #pylab.show()

    if useLabels and labels.L is not None :
	numPatterns = 0
	for i in range(labels.numClasses) :
	    numPatterns += labels.classSize[i]
	    #pylab.figtext(0.05, float(numPatterns) / len(labels), labels.classLabels[i])
	    #pylab.figtext(float(numPatterns) / len(labels), 0.05, labels.classLabels[i])
	    pylab.axhline(numPatterns, color = 'black', linewidth = 1)
	    pylab.axvline(numPatterns, color = 'black', linewidth = 1)
    pylab.axis([0, len(labels), 0, len(labels)])
    if fileName is not None :
        pylab.savefig(fileName)
	pylab.close()

def my_lines(ax, pos, *args, **kwargs):
    if ax == 'x':
        for p in pos:
            plt.axvline(p, *args, **kwargs)
    else:
        for p in pos:
            plt.axhline(p, *args, **kwargs)

def transform(gameData):
# refHistory= { ref: { team : [total, games] }
    winData = transformWinLoss(gameData)
    collapsedData = winData.groupby([winData.referee,winData.team]).sum()
    collapsedData['winPercentage'] = collapsedData["teamWin"]/collapsedData["teamPlay"]
    totalData = winData.groupby(winData.team).sum()
    totalData['winPercentage'] = totalData["teamWin"]/totalData["teamPlay"]
    print totalData.sort_index(by='winPercentage', ascending=False)
    plt.figure(figsize=(8, 8))
    #totalData['winPercentage'].plot(kind='bar', alpha=0.5)
    #plt.show()
    #plt.clf()
    print totalData.index.values
    for team in totalData.index.values:
        print "doing %s" % team
        teamData = winData.loc[winData['team'] == team ].groupby(winData.referee).sum()
        teamData = teamData[teamData['teamPlay'] > 6]
        if teamData.empty:
            continue
        print teamData
        teamData['winPercentage'] = teamData["teamWin"]/teamData["teamPlay"]
        teamData['winPercentage'].plot(kind='bar', alpha=0.5)
        teamAvg = totalData['winPercentage'][team]
        print teamAvg
        plt.ylim([0,1.0])
        plt.axhline(teamAvg, color='k')
        plt.savefig('%s.png' % team)
        plt.clf()

def mark_cross(center, **kwargs):
    """Mark a cross. Correct for matplotlib imshow funny coordinate system.
    """
    N = 20
    plt.hold(1)
    plt.axhline(y=center[1]-0.5, **kwargs)
    plt.axvline(x=center[0]-0.5, **kwargs)

 def plot_eigen_spectrum(eig, laser_freq_eV=0.0, file_to_save='eig.png', figuresize=(10,10)):
     plt.figure(figsize=figuresize) # in inches!
     plt.plot(eig.real, 'o', ms=2)
     plt.ylabel(r'$E, eV$', fontsize=26)
     plt.axhline(y=laser_freq_eV, color='r', ls='--')
     plt.axhline(y=0, color='r', ls='--')
     if file_to_save:
         plt.savefig(file_to_save)
         plt.close()
     return None

def epi_vs_gain_volcano_plot(filtered_gain_snps, filtered_epi_snps, gain_vals, epi_vals, max_p, min_I3):
    gain_I3 = []
    gain_log_p = []
    for snps in filtered_gain_snps:
        gain_I3.append(gain_vals[snps])

        order = switch_snp_key_order(snps)
        if epi_vals.has_key(order[0]):
            gain_log_p.append(epi_vals[order[0]])
        elif epi_vals.has_key(order[1]):
            gain_log_p.append(epi_vals[order[1]])
    gain_log_p = -1 * np.log10(gain_log_p)

    epi_I3 = []
    epi_log_p = []
    for snps in filtered_epi_snps:
        order = switch_snp_key_order(snps)
        if gain_vals.has_key(order[0]):
            epi_I3.append(gain_vals[order[0]])
        elif gain_vals.has_key(order[1]):
            epi_I3.append(gain_vals[order[1]])

        epi_log_p.append(epi_vals[snps])
    epi_log_p = -1 * np.log10(epi_log_p)

    mp.figure(1)
    mp.xlabel("I3")
    mp.ylabel("-log10(P)")
    mp.title("Volcano plot - EPISTASIS and GAIN")
    mp.plot(epi_I3, epi_log_p, "bo")
    mp.plot(gain_I3, gain_log_p, "ro")
    mp.axhline(y=(-1 * np.log10(max_p)), linewidth=2, color="g")
    mp.axvline(x=min_I3, linewidth=2, color="g")
    # label max point
    max_x = np.max(gain_I3)
    max_y = np.max(gain_log_p)
    best_connection = ""
    # label best edge
    for snps in epi_vals:
        if -1 * np.log10(epi_vals[snps]) == max_y:
            best_connection = str(snps)
    mp.text(
        np.max(gain_I3),
        np.max(gain_log_p),
        best_connection,
        fontsize=10,
        horizontalalignment="center",
        verticalalignment="center",
    )
    mp.show()

    print

def plot_episode(args):
    """Plot an episode plucked from the large h5 database"""
    print "plot_episode"
    # load the data file
    tblfilename = "bf_optimize_mavlink.h5"
    h5file = tb.open_file(tblfilename, mode = "a")
    # get the table handle
    table = h5file.root.v2.evaluations

    # selected episode
    episode_row = table.read_coordinates([int(args.epinum)])
    # compare episodes
    episode_row_1 = table.read_coordinates([2, 3, 22, 46]) # bad episodes
    print "row_1", episode_row_1.shape
    # episode_row = table.read_coordinates([3, 87])
    episode_target = episode_row["alt_target"]
    episode_target_1 = [row["alt_target"] for row in episode_row_1]
    print "episode_target_1.shape", episode_target_1
    episode_timeseries = episode_row["timeseries"][0]
    episode_timeseries_1 = [row["timeseries"] for row in episode_row_1]
    print "row", episode_timeseries.shape
    print "row_1", episode_timeseries_1

    sl_start = 0
    sl_end = 2500
    sl_len = sl_end - sl_start
    sl = slice(sl_start, sl_end)
    pl.plot(episode_timeseries[sl,1], "k-", label="alt", lw=2.)
    print np.array(episode_timeseries_1)[:,:,1]
    pl.plot(np.array(episode_timeseries_1)[:,:,1].T, "k-", alpha=0.2)
    # alt_hold = episode_timeseries[:,0] > 4
    alt_hold_act = np.where(episode_timeseries[sl,0] == 11)
    print "alt_hold_act", alt_hold_act[0].shape, sl_len
    alt_hold_act_min = np.min(alt_hold_act)
    alt_hold_act_max = np.max(alt_hold_act)
    print "min, max", alt_hold_act_min, alt_hold_act_max, alt_hold_act_min/float(sl_len), alt_hold_act_max/float(sl_len),

    # pl.plot(episode_timeseries[sl,0] * 10, label="mode")
    pl.axhspan(-100., 1000,
               alt_hold_act_min/float(sl_len),
               alt_hold_act_max/float(sl_len),
               facecolor='0.5', alpha=0.25)
    pl.axhline(episode_target, label="target")
    pl.xlim((0, sl_len))
    pl.xlabel("Time steps [1/50 s]")
    pl.ylabel("Alt [cm]")
    pl.legend()
    if args.plotsave:
        pl.gcf().set_size_inches((10, 3))
        pl.gcf().savefig("%s.pdf" % (sys.argv[0][:-3]), dpi=300, bbox_inches="tight")
    pl.show()

def diff_plot_bar(lists, list_ids, xticks,
                  rotation=0, xlabel='', ylabel='Accuracy',
                  hline_at=None, legend_title='Method', **kwargs):
    """
    Compare the scores of paired of experiment ids and plot a bar chart or their accuracies.
    :param list1, list2: [1,2,3], [4,5,6] means exp 1 is compared to exp 4, etc ...
    :param list1_id, list2_id: name for the first/ second group of experiments, will appear in legend
    :param xticks: labels for the x-axis, one per pair of experiments, e.g.
    list('abc') will label the first pair 'a', etc. Will appear as ticks on x axis.
    If only two lists are provided a significance test is run for each pair and a * is added if pair is
    significantly different
    :param rotation: angle of x axis ticks
    :param hline_at: draw a horizontal line at y=hline_at. Useful for baselines, etc
    :param kwargs: extra arguments for sns.factorplot
    """
    assert len(set(map(len, lists))) == 1
    assert len(list_ids) == len(lists)

    df_scores, df_reps, df_groups, df_labels = [], [], [], []
    if len(lists) == 2:
        for i, (a, b) in enumerate(zip(*lists)):
            significance_df, names, mean_scores = get_demsar_params([a, b],
                                                                    name_format=['id',
                                                                                 'expansions__vectors__id',
                                                                                 'expansions__vectors__composer',
                                                                                 'expansions__vectors__algorithm',
                                                                                 'expansions__vectors__dimensionality'])
            if significance_df is None:
                continue
            if significance_df.significant[0]:
                xticks[i] += '*'

    for i, exp_ids in enumerate(zip(*lists)):
        data, folds = get_cv_scores_many_experiment(exp_ids)
        df_scores.extend(data)
        df_reps.extend(folds)
        df_labels.extend(len(folds) * [xticks[i]])
        for list_id in list_ids:
            df_groups.extend(len(folds) // len(lists) * [list_id])

    df = pd.DataFrame(dict(Accuracy=df_scores, reps=df_reps, Method=df_groups, labels=df_labels))
    df.rename(columns={'Method': legend_title}, inplace=True)
    g = sns.factorplot(y='Accuracy', hue=legend_title, x='labels', data=df, kind='bar', aspect=1.5, **kwargs);
    g.set_xticklabels(rotation=rotation);
    # remove axis labels
    for ax in g.axes.flat:
        ax.set(xlabel=xlabel, ylabel=ylabel)
    if hline_at is not None:
        plt.axhline(hline_at, color='black')

def plot_moving_correlation(series, chrons, start, end, 
                            crit=None, window=get_window('boxcar', 31)):

  w = len(window) / 2
  t = range(start + w, end - w + 1)
  
  p = series.ix[start:end].values.flatten()
 
  for col in chrons:
    c = chrons[col].ix[start:end].values
    corrs, pvals = moving_correlation(p, c, window)
    plt.plot(t, corrs, label=col.upper(), markevery=(t[-1]-t[0])/12, **pens[col.upper()])

  if crit is not None:
    plt.axhline(y=crit, linestyle='--', color='black')  

def plot_all_pair_correlation(dirname=".",target="g_tilde"):
    
    handles=[]
    global folders
    plt.figure(figsize=(20,10))
    global folders
    # plot all energies
    print "Pair correlations"
    for f in folders:
        label=get_label(target,f)
        handles.append(anal.plot_file(f.dir_path+"/pair_correlation_t.dat",label=label,error=True ))
        
    plt.legend(handles=handles)
    plt.axhline(y=1, xmin=0,hold=True)
    
    handles=[]

def plot_weather(p, ranges, plot_happiness = False, threshold=40.0):

    plts = []

    for date_start, date_end in ranges:
        print "Getting weather data for %s - %s..." % (date_start, date_end)
        station_info, observations = p.get_weather_hourly(94846, date_start, date_end)
        observations = Plenario.get_weather_observations_list(observations, "datetime", ["drybulb_fahrenheit"])
        observations.reverse()

        dates = [x[0] for x in observations]
        datesnorm = normalize_dates(dates)
        temps = [x[1] for x in observations]
        fill_missing_values(temps)


        label = "%s - %s" % (date_start, date_end)

        if plot_happiness:
            Ys = gen_happiness(temps, threshold)
        else:
            Ys = temps

        w_plt, = plt.plot(datesnorm, Ys, label=label)
        plts.append(w_plt)

    if plot_happiness:
        hline = 0
        title = "Weather-related Happiness in Chicago"
        ylabel = "Happiness"
    else:
        hline = 32
        title = "Temperatures in Chicago"
        ylabel = u"Temperature (°F)"

    xlabel = "Week"

    nweeks = (int(datesnorm[-1]) / 604800) + 1

    plt.xticks([w*604800 for w in range(nweeks)], [w for w in range(nweeks)])

    plt.axhline(hline, color="gray", linestyle="--")
    plt.legend(handles=plts, loc="lower left")
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.show()

def plotDifference(bst, dtest, target = None):
    p = bst.predict(dtest)
    l = 0
    if (target is None):    
        l = dtest.get_label()
    else:
        l = target
    res = np.abs((l-p)/l)
    res = res[~np.isnan(res)]
    res = res[~np.isinf(res)]
    xmax = 0.2
    for y in np.arange(0,xmax,xmax/8):    
        plt.axvline(y, color='g', linestyle='dashed')
    plt.hist(res, range=(0,xmax), bins=50, color='red')
    plt.axvline(res.mean(), color='c', linestyle='dashed', linewidth=2)
    for y in range(200,1600,200):    
        plt.axhline(y, color='g', linestyle='dashed')

 def draw_cost_history(self, filename='output/cost.png', compare_dict=None):
     colors = ['blue', 'red', 'yellow', 'green', 'pink']
     f, ax = plt.subplots()
     ax.plot(self.cost_history, lw=1, c='black', label='cost', alpha=0.8)
     df = pd.DataFrame(columns=['rolling mean'], data=self.cost_history)
     df['rolling mean'] = pd.rolling_mean(df['rolling mean'], window=100)
     df.plot(ax=ax, lw=3, c='green')
     if compare_dict is not None:
         for idx, (name, vals) in enumerate(compare_dict.iteritems()):
             c = colors[idx % len(colors)]
             if hasattr(vals, '__iter__'):
                 plt.plot(vals, lw=2, label=name, c=c, alpha=0.6)
             else:
                 plt.axhline(y=vals, lw=2, label=name, c=c, alpha=0.6)
     plt.legend(loc='best')
     plt.savefig(filename, dpi=150)
     plt.close('all')

def view_design(g, t, arange):
    """Plot the scaling and wavelet kernel.

    Plot the input scaling function and wavelet kernels, indicates the wavelet
    scales by legend, and also shows the sum of squares G and corresponding
    frame bounds for the transform.

    Parameters
    ----------
    g : list of  function handles for scaling function and wavelet kernels
    t : array of wavelet scales corresponding to wavelet kernels in g
    arange : approximation range

    Returns
    -------
    h : figure handle
    """
    x = np.linspace(arange[0], arange[1], 1e3)
    h = plt.figure()
    
    J = len(g) 
    G = np.zeros(x.size)

    for n in range(J):
        if n == 0:
            lab = 'h'
        else:
            lab = 't=%.2f' % t[n-1]
        plt.plot(x, g[n](x), label=lab)
        G += g[n](x)**2

    plt.plot(x, G, 'k', label='G')

    (A, B, _, _) = framebounds(g, arange[0], arange[1])
    plt.axhline(A, c='m', ls=':', label='A')
    plt.axhline(B, c='g', ls=':', label='B')
    plt.xlim(arange[0], arange[1])

    plt.title('Scaling function kernel h(x), Wavelet kernels g(t_j x) \n'
              'sum of Squares G, and Frame Bounds')
    plt.yticks(np.r_[0:4])
    plt.ylim(0, 3)
    plt.legend()

    return h

def simSeries(numSeries):
    prob = 0.5
    fracWon = []
    probs = []
    while prob <= 1.0:
        seriesWon = 0
        for i in range(numSeries):
            if playSeries(7, prob):
                seriesWon += 1
        fracWon.append(seriesWon/numSeries)
        probs.append(prob)
        prob += 0.01
    plt.plot(probs, fracWon, linewidth=3)
    plt.xlabel("Probability of Winning a Game")
    plt.ylabel("Probability of Winning a Series")
    plt.axhline(0.95)
    plt.ylim(0.5, 1.1)
    plt.title(str(numSeries) + " Seven-Game Series")