def plot_degreeRate(db, keynames, save_path):
	degRate_x_name = 'degRateDistr_x'
	degRate_y_name = 'degRateDistr_y'

	plt.clf()
	plt.figure(figsize = (8, 5))

	plt.subplot(1, 2, 1)
	plt.plot(db[keynames['mog']][degRate_x_name], db[keynames['mog']][degRate_y_name], 'b-', lw=5, label = 'fairyland')
	plt.plot(db[keynames['mblg']][degRate_x_name], db[keynames['mblg']][degRate_y_name], 'r:', lw=5, label = 'twitter')
	plt.plot(db[keynames['im']][degRate_x_name], db[keynames['im']][degRate_y_name], 'k--', lw=5, label = 'yahoo')
	plt.xscale('log')
	plt.grid(True)
	plt.title('interaction')
	plt.legend(('fairyland', 'twitter', 'yahoo'), loc = 4, prop = {'size': 10})
	plt.xlabel('In-degree to Out-degree Ratio')
	plt.ylabel('CDF')

	plt.subplot(1, 2, 2)
	plt.plot(db[keynames['mogF']][degRate_x_name], db[keynames['mogF']][degRate_y_name], 'b-', lw=5, label = 'fairyland')
	plt.plot(db[keynames['mblgF']][degRate_x_name], db[keynames['mblgF']][degRate_y_name], 'r:', lw=5, label = 'twitter')
	#plt.plot(db[keynames['imF']][degRate_x_name], db[keynames['imF']][degRate_y_name], 'k--', lw=5, label = 'yahoo')
	plt.xscale('log')
	plt.grid(True)
	plt.title('ally')
	plt.xlabel('In-degree to Out-degree Ratio')
	plt.ylabel('CDF')

	plt.savefig(os.path.join(save_dir, save_path))

def plot_feat_hist(data_name_list, filename=None):
    pylab.clf()
    # import pdb;pdb.set_trace()
    num_rows = 1 + (len(data_name_list) - 1) / 2
    num_cols = 1 if len(data_name_list) == 1 else 2
    pylab.figure(figsize=(5 * num_cols, 4 * num_rows))

    for i in range(num_rows):
        for j in range(num_cols):
            pylab.subplot(num_rows, num_cols, 1 + i * num_cols + j)
            x, name = data_name_list[i * num_cols + j]
            pylab.title(name)
            pylab.xlabel('Value')
            pylab.ylabel('Density')
            # the histogram of the data
            max_val = np.max(x)
            if max_val <= 1.0:
                bins = 50
            elif max_val > 50:
                bins = 50
            else:
                bins = max_val
            n, bins, patches = pylab.hist(
                x, bins=bins, normed=1, facecolor='green', alpha=0.75)

            pylab.grid(True)

    if not filename:
        filename = "feat_hist_%s.png" % name

    pylab.savefig(os.path.join(CHART_DIR, filename), bbox_inches="tight")

    def __init__(self, fig, gs, num_plots, rows=None, cols=None):
        if cols is None:
            cols = int(np.floor(np.sqrt(num_plots)))
        if rows is None:
            rows = int(np.ceil(float(num_plots)/cols))
        assert num_plots <= rows*cols, 'Too many plots to put into gridspec.'

        self._fig = fig
        self._gs = gridspec.GridSpecFromSubplotSpec(8, 1, subplot_spec=gs)
        self._gs_legend = self._gs[0:1, 0]
        self._gs_plot   = self._gs[1:8, 0]

        self._ax_legend = plt.subplot(self._gs_legend)
        self._ax_legend.get_xaxis().set_visible(False)
        self._ax_legend.get_yaxis().set_visible(False)

        self._gs_plots = gridspec.GridSpecFromSubplotSpec(rows, cols, subplot_spec=self._gs_plot)
        self._axarr = [plt.subplot(self._gs_plots[i], projection='3d') for i in range(num_plots)]
        self._lims = [None for i in range(num_plots)]
        self._plots = [[] for i in range(num_plots)]

        for ax in self._axarr:
            ax.tick_params(pad=0)
            ax.locator_params(nbins=5)
            for item in (ax.get_xticklabels() + ax.get_yticklabels() + ax.get_zticklabels()):
                item.set_fontsize(10)

        self._fig.canvas.draw()
        self._fig.canvas.flush_events()   # Fixes bug with Qt4Agg backend

    def plot_histogram(self, main="", numrows=1, numcols=1, fignum=1):
        """Plot a histogram of choices and probability sums. Expects probabilities as (at least) a 2D array.
        """
        from matplotlib.pylab import bar, xticks, yticks, title, text, axis, figure, subplot

        probabilities = self.get_probabilities()
        if probabilities.ndim < 2:
            raise StandardError, "probabilities must have at least 2 dimensions."
        alts = probabilities.shape[1]
        width_par = (1 / alts + 1) / 2.0
        choice_counts = self.get_choice_histogram(0, alts)
        sum_probs = self.get_probabilities_sum()

        subplot(numrows, numcols, fignum)
        bar(arange(alts), choice_counts, width=width_par)
        bar(arange(alts) + width_par, sum_probs, width=width_par, color="g")
        xticks(arange(alts))
        title(main)
        Axis = axis()
        text(
            alts + 0.5,
            -0.1,
            "\nchoices histogram (blue),\nprobabilities sum (green)",
            horizontalalignment="right",
            verticalalignment="top",
        )

    def plot_corner_posteriors(self, savefile=None, labels=["T1", "R1", "Av", "T2", "R2"]):
        '''
        Plots the corner plot of the MCMC results.
        '''
        ndim = len(self.sampler.flatchain[0,:])
        chain = self.sampler
        samples = chain.flatchain
        
        samples = samples[:,0:ndim]  
        plt.figure(figsize=(8,8))
        fig = corner.corner(samples, labels=labels[0:ndim])
        plt.title("MJD: %.2f"%self.mjd)
        name = self._get_save_path(savefile, "mcmc_posteriors")
        plt.savefig(name)
        plt.close("all")
        

        plt.figure(figsize=(8,ndim*3))
        for n in range(ndim):
            plt.subplot(ndim,1,n+1)
            chain = self.sampler.chain[:,:,n]
            nwalk, nit = chain.shape
            
            for i in np.arange(nwalk):
                plt.plot(chain[i], lw=0.1)
                plt.ylabel(labels[n])
                plt.xlabel("Iteration")
        name_walkers = self._get_save_path(savefile, "mcmc_walkers")
        plt.tight_layout()
        plt.savefig(name_walkers)
        plt.close("all")  

def EnhanceContrast(g, r=3, op_kernel=15, silence=True):
    
    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(op_kernel,op_kernel))
    opening = cv2.morphologyEx(g, cv2.MORPH_OPEN, kernel)
    
    g_copy = np.asarray(np.copy(g), dtype=np.float)

    m_f = np.mean(opening)
        
    u_max = 245; u_min = 10; t_min = np.min(g); t_max = np.max(g)

    idx_gt_mf = np.where(g_copy > m_f)
    idx_lt_mf = np.where(g_copy <= m_f)

    g_copy[idx_gt_mf] = -0.5 * ((u_max-u_min) / (m_f-t_max)**r) * (g_copy[idx_gt_mf]-t_max)**r + u_max
    g_copy[idx_lt_mf] = 0.5 * ((u_max-u_min) / (m_f-t_min)**r) * (g_copy[idx_lt_mf]-t_min)**r + u_min 

    if silence == False:
        plt.subplot(1,2,1)
        plt.imshow(g, cmap='gray')
        plt.title('Original image')
        plt.subplot(1,2,2)
        plt.imshow(g_copy, cmap='gray')
        plt.title('Enhanced image')
        plt.show()
        
    return g_copy

def test_likelihood_evaluator3():
    
    tr = template.TemplateRenderCircleBorder()
    tr.set_params(14, 6, 4)

    t1 = tr.render(0, np.pi/2)
    img = np.zeros((240, 320), dtype=np.uint8)

    env = util.Environmentz((1.5, 2.0), (240, 320))
    
    le2 = likelihood.LikelihoodEvaluator3(env, tr)

    img[(120-t1.shape[0]/2):(120+t1.shape[0]/2), 
        (160-t1.shape[1]/2):(160+t1.shape[1]/2)] += t1 *255
    pylab.subplot(1, 2, 1)
    pylab.imshow(img, interpolation='nearest', cmap=pylab.cm.gray)

    state = np.zeros(1, dtype=util.DTYPE_STATE)

    xvals = np.linspace(0, 2.,  100)
    yvals = np.linspace(0, 1.5, 100)
    res = np.zeros((len(yvals), len(xvals)), dtype=np.float32)
    for yi, y in enumerate(yvals):
        for xi, x in enumerate(xvals):
            state[0]['x'] = x
            state[0]['y'] = y
            state[0]['theta'] = np.pi / 2. 
            res[yi, xi] =     le2.score_state(state, img)
    pylab.subplot(1, 2, 2)
    pylab.imshow(res)
    pylab.colorbar()
    pylab.show()

 def plot(x,y,field,filename,c=200):
     plt.figure()
     # define grid.
     xi = np.linspace(min(x),max(x),100)
     yi = np.linspace(min(y),max(y),100)
     # grid the data.
     si_lin = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     si_cub = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     print np.min(field)
     print np.max(field)
     plt.subplot(211)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_lin,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_lin,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Lineaarinen interpolointi')
     #plt.tight_layout()
     plt.subplot(212)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_cub,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_cub,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Kuubinen interpolointi')
     plt.savefig(filename)

def fit_calib_auto(X_bh, X_bg_gh, do_l1=False):
    """
Does auto-regression on the 6-DOF variables : x,y,z,r,p,y.
"""
    assert X_bh.shape[0]==X_bg_gh.shape[0]==12, "calib data has unknown shape."
    
    #X_bh = X_bh[:,500:1000]
    #X_bg_gh = X_bg_gh[:,500:1000]
    
    axlabels = ['x','y','z','roll','pitch','yaw']
    for k in xrange(1,100, 10):
        print "order : k=", k
        plt.clf()
        W = np.empty((6, 2*k+1))
        for i in xrange(6):
            j = i+3 if i > 2 else i
            Ai, bi, W[i,:] = fit_auto(X_bh[j,:], X_bg_gh[j,:], k, do_l1)
            est = Ai.dot(W[i,:])
            print " norm err : ", np.linalg.norm(bi - est)
    
            plt.subplot(3,2,i+1)
            plt.plot(bi, label='pr2')
            plt.plot(X_bh[j,k-1:], label='hydra')
            plt.plot(est, label='estimate')
            plt.ylabel(axlabels[i])
            plt.legend()
        plt.show()

def PlotMtxError(Corr_w):
    max_val = 1
    min_val = -0.1

    AvCorr = np.sum(Corr_w, axis=0)
    dCorr = Corr_w - AvCorr
    errCorr = np.log10(np.sqrt(np.einsum("i...,i...", dCorr, dCorr)) / np.absolute(AvCorr) / np.sqrt(Corr_w.shape[0]))
    # print errCorr.shape
    # print errCorr

    plt.rcParams.update({"font.size": 6, "font.weight": "bold"})
    for i in xrange(errCorr.shape[0]):
        plt.subplot(2, 7, i + 1)
        plt.title("SITE " + str(i + 1) + ":: \nHistogram of errors in corr. mtx.")
        plt.hist(errCorr[0, :, :].flatten(), 256, range=(min_val, max_val))
        plt.xlabel("log_10(sigma)")
        plt.ylabel("Count")

        plt.subplot(2, 7, i + 7 + 1)
        plt.imshow(errCorr[0, :, :], vmin=min_val, vmax=max_val)
        cbar = plt.colorbar(shrink=0.25, aspect=40)
        cbar.set_label("log_10(sigma)")
        plt.set_cmap("gist_yarg")
        plt.title("SITE " + str(i + 1) + ":: \nError in corr. matx. values")
        plt.xlabel("Site i")
        plt.ylabel("Site j")
    plt.show()

def plot_hyperplane(X, Y, model, K, plot_id, d = 500):
    I0 = np.where(Y==-1)[0]
    I1 = np.where(Y==1)[0]

    plt.subplot(plot_id)

    plt.plot(X[I1, 0], X[I1, 1], 'og')
    plt.plot(X[I0, 0], X[I0, 1], 'xb')

    min_val = np.min(X, 0)
    max_val = np.max(X, 0)

    clf = model()
    clf.train(X, Y, K)

    x0_plot = np.linspace(min_val[0, 0], max_val[0, 0], d)
    x1_plot = np.linspace(min_val[0, 1], max_val[0, 1], d)

    [x0, x1] = plt.meshgrid(x0_plot, x1_plot);

    Y_all = np.matrix(np.zeros([d, d]))

    for i in range(d):
        X_all = np.matrix(np.zeros([d, 2]))
        X_all[:, 0] = np.matrix(x0[:, i]).T
        X_all[:, 1] = np.matrix(x1[:, i]).T
        Y_all[:, i] = clf.predict(X_all)

    plt.contour(np.array(x0), np.array(x1), np.array(Y_all), levels = [0.0], colors = 'red')

def plot_reconstruction_result(res):
    """ Plot original and reconstructed graph plus time series
    """
    fig = plt.figure(figsize=(32, 8))
    gs = mpl.gridspec.GridSpec(1, 4)

    # original graph
    orig_ax = plt.subplot(gs[0])
    plot_graph(nx.from_numpy_matrix(res.A.orig), orig_ax)
    orig_ax.set_title('Original graph')

    # time series
    ax = plt.subplot(gs[1:3])
    sns.tsplot(
        time='time', value='theta',
        unit='source', condition='oscillator',
        estimator=np.mean, legend=False,
        data=compute_solutions(res),
        ax=ax)
    ax.set_title(r'$A_{{err}} = {:.2}, B_{{err}} = {:.2}$'.format(*compute_error(res)))

    # reconstructed graph
    rec_ax = plt.subplot(gs[3])
    tmp = res.A.rec
    tmp[abs(tmp) < 1e-1] = 0
    plot_graph(nx.from_numpy_matrix(tmp), rec_ax)
    rec_ax.set_title('Reconstructed graph')

    plt.tight_layout()
    save(fig, 'reconstruction_overview')

    def train(self, ratings, model_path):
 
        self.mu = ratings.mean
        self.P = 0.001 * np.matrix(np.random.randn(len(ratings.rows), self.num_factor))
        self.bu = 0.001 * np.matrix(np.random.randn(len(ratings.rows), 1))
        self.Q = 0.001 * np.matrix(np.random.randn(len(ratings.cols), self.num_factor))
        self.bi = 0.001 * np.matrix(np.random.randn(len(ratings.cols), 1))
        
        self.rows = dict(ratings.rows)
        self.cols = dict(ratings.cols)       

        if self.validate > 0:
            T = ratings.kv_dict.items()
            random.shuffle(T)
            k = len(T) / self.validate
            self.L_validate = T[0 : k]
            self.L_train = T[k :]
        else:
            self.L_train = ratings.kv_dict.items()

        rmse_train = [0.0] * self.max_iter
        rmse_validate = [0.0] * self.max_iter       
 
        for s in range(self.max_iter):

            random.shuffle(self.L_train)
            self.current_sample = 0
            self.sqr_err = 0.0

            self.threads = [ParallelSGD('Thread_%d' % n, self) for n in range(self.num_thread)]
            
            start = time.time()
            for t in self.threads:
                t.start()
                t.join()
            terminal = time.time()

            duration = terminal - start

            rmse_train[s] = math.sqrt(self.sqr_err / len(ratings.kv_dict))
    
            if self.validate > 0:
                m = SparseMatrix()
                m.kv_dict = {k : v for (k, v) in self.L_validate}
                rmse_validate[s] = float(self.test(m))
            
            sys.stderr.write('Iter: %4.4i' % (s + 1))
            sys.stderr.write('\t[Train RMSE] = %f' % rmse_train[s])
            if self.validate > 0:
                sys.stderr.write('\t[Validate RMSE] = %f' % rmse_validate[s])
            sys.stderr.write('\t[Duration] = %f' % duration)
            sys.stderr.write('\t[Samples] = %d\n' % len(self.L_train))

            self.dump_model(model_path + '/' + 'model_%4.4i' % (s + 1))
            self.dump_raw_model(model_path + '/' + 'model_%4.4i.raw_model' % (s + 1))

        plt.subplot(111)
        plt.plot(range(self.max_iter), rmse_train, '-og')
        plt.plot(range(self.max_iter), rmse_validate, '-xb')
        plt.show()

 def sanity_example2(self):
   """
     Checking Kepler orbit calculation example.
   """
   import numpy
   from PyAstronomy import pyasl
   import matplotlib.pylab as plt
   
   # Instantiate a Keplerian elliptical orbit with
   # semi-major axis of 1.3 length units,
   # period of 2 time units, eccentricity of 0.5, and
   # longitude of ascending node of 70 degrees.
   ke = pyasl.KeplerEllipse(1.3, 2., e=0.5, Omega=70.)
   
   # Get a time axis
   t = numpy.linspace(0, 6.5, 200)
   
   # Calculate the orbit position at the given points
   # in a Cartesian coordinate system.
   pos = ke.xyzPos(t)
   print "Shape of output array: ", pos.shape
   
   # x, y, and z coordinates for 50th time point
   print "x, y, z for 50th point: ", pos[50, ::]
   
   # Calculate orbit radius as a function of the
   radius = ke.radius(t)
   
   # Plot x and y coordinates of the orbit
   plt.subplot(2,1,1)
   plt.plot(pos[::,0], pos[::,1], 'bp')
   # Plot orbit radius as a function of time
   plt.subplot(2,1,2)
   plt.plot(t, radius, 'bp')

def test():
	## Load files
    s = load_spectrum('ring28yael')
    w = linspace(1510e-9,1600e-9,len(s))
    
	## Process
    mins = find_minima(s)
    w_p = 1510e-9 + array(mins) * 90.e-9/len(w)
    ww = 2 * pi * 3e8/w_p   
    
	## Plot
    pl.plot(w,s)
    pl.plot(w_p,s[mins],'o')
    pl.show()
    
    beta2 = -1./(112e-6*2*pi)*diff(diff(ww))/(diff(ww)[:-1]**3)
    p = polyfit(w_p[1:-1], beta2, 1)
    
    savetxt('ring28yael-p.txt', w_p)
    
    pl.subplot(211)
    pl.plot(w,s)
    pl.plot(w_p,s[mins],'o')
    
    pl.subplot(212)
    pl.plot(w_p[1:-1]*1e6, beta2)
    pl.plot(w_p[1:-1]*1e6, p[1]+ p[0]*w_p[1:-1], label="q=%.2e"%p[0])
    pl.legend()
        
    pl.show()

def chooseDegree(npts, mindegree=0, maxdegree=20, filename=None):
    """Gets noisy data, uses cross validation to estimate error, and fits new data with best model."""
    x, y = bv.noisyData(npts)
    degrees = numpy.arange(mindegree,maxdegree+1)
    errs = numpy.zeros_like(degrees,dtype=numpy.float)
    for i,d in enumerate(degrees):
        errs[i] = estimateError(x, y, d)

    plt.subplot(1,2,1)
    plt.plot(degrees,errs,'bo-')
    plt.xlabel("Degree")
    plt.ylabel("CV Error")

    besti = numpy.argmin(errs)
    bestdegree = degrees[besti]

    plt.subplot(1,2,2)
    x2, y2 = bv.noisyData(npts)
    plt.plot(x2,y2,'ro')
    xs = numpy.linspace(min(x),max(x),150)
    fitf = numpy.poly1d(numpy.polyfit(x2,y2,bestdegree))
    plt.plot(xs,fitf(xs),'g-',lw=2)
    plt.xlim((bv.MIN,bv.MAX))
    plt.ylim((-2.,2.))
    plt.suptitle('Selected Degree '+str(bestdegree))
    bv.outputPlot(filename)

def find_params():

    FRAMES =  np.arange(30)*100

    frame_images = organizedata.get_frames(ddir("bukowski_04.W2"), FRAMES)
    print "DONE READING DATA"

    CLUST_EPS = np.linspace(0, 0.5, 10)
    MIN_SAMPLES = [2, 3, 4, 5]
    MIN_DISTS = [2, 3, 4, 5, 6]
    THOLD = 240

    fracs_2 = np.zeros((len(CLUST_EPS), len(MIN_SAMPLES), len(MIN_DISTS)))

    for cei, CLUST_EP in enumerate(CLUST_EPS):
        for msi, MIN_SAMPLE in enumerate(MIN_SAMPLES):
            for mdi, MIN_DIST in enumerate(MIN_DISTS):
                print cei, msi, mdi
                numclusters = np.zeros(len(FRAMES))
                for fi, im in enumerate(frame_images):
                    centers = frame_clust_points(im, THOLD, MIN_DIST, 
                                                 CLUST_EP, MIN_SAMPLE)
                    # cluster centers
                    numclusters[fi] = len(centers)
                fracs_2[cei, msi, mdi] = float(np.sum(numclusters == 2))/len(numclusters)
    pylab.figure(figsize=(12, 8))
    for mdi, MIN_DIST in enumerate(MIN_DISTS):
        pylab.subplot(len(MIN_DISTS), 1, mdi+1)

        for msi in range(len(MIN_SAMPLES)):
            pylab.plot(CLUST_EPS, fracs_2[:, msi, mdi], label='%d' % MIN_SAMPLES[msi])
        pylab.title("min_dist= %3.2f" % MIN_DIST)
    pylab.legend()
    pylab.savefig('test.png', dpi=300)

def plot_env(env, growth_data=None, log_pop_size=True):
    """
    Plot environment and the population growth (if given).
    """
    if growth_data is None:
        df = env.as_df(melted=True)
        g = sns.factorplot(x="index", y="value", hue="nutrient", data=df)
        return g
    # get un-melted dataframe
    df = env.as_df(melted=False)
    df["t"] = growth_data["t"].t
    plt.subplot(2, 1, 1)
    df["pop_size"] = np.exp(growth_data["log_pop_size"])
    for nutr in env.nutrs:
        plt.plot(df["t"], df[nutr], label=nutr)
    plt.legend()
    plt.xlabel("Time")
    plt.subplot(2, 1, 2)
    plt.xlabel("Time")
    if log_pop_size:
        plt.plot(df["t"], np.log2(df["pop_size"]))
        plt.ylabel("Pop. size (log$_\mathrm{2}$)")
    else:
        plt.plot(df["t"], df["pop_size"])
        plt.ylabel("Pop. size")

 def ACF_PACF_plot(self):
     #plot ACF and PACF to find the number of terms needed for the AR and MA in ARIMA
     # ACF finds MA(q): cut off after x lags 
     # and PACF finds AR (p): cut off after y lags 
     # in ARIMA(p,d,q) 
     lag_acf = acf(self.ts_log_diff, nlags=20)
     lag_pacf = pacf(self.ts_log_diff, nlags=20, method='ols')
     
     #Plot ACF:
     ax=plt.subplot(121)
     plt.plot(lag_acf)
     ax.set_xlim([0,5])
     plt.axhline(y=0,linestyle='--',color='gray')
     plt.axhline(y= -1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
     plt.axhline(y= 1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
     plt.title('Autocorrelation Function')
     
     #Plot PACF:
     plt.subplot(122)
     plt.plot(lag_pacf)
     plt.axhline(y=0,linestyle='--',color='gray')
     plt.axhline(y= -1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
     plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
     plt.title('Partial Autocorrelation Function')
     plt.tight_layout()

def plot_all_dep(data, with_subplots=True, remove_weekday=False):
    country_codes = list(data.country_code.unique())
    channels = list(data.marketing_channel.unique())
    ny = len(country_codes)
    nx = len(channels)
    plot_num = 1

    abrevs = channel_abbrevs()

    for channel in channels:
        for country_code in country_codes:
            data_one = data[(data.country_code == country_code)
                            & (data.marketing_channel == channel)]
            if len(data_one) > 2:
                if with_subplots:
                    plt.subplot(nx, ny, plot_num)
                dates = data_one.date
                y = data_one.user_visits
                if remove_weekday:
                    y = remove_weekday_seasonality(dates, y)
                plt.plot(dates, y)
                if with_subplots:
                    title = abrevs[channel] + '_' + country_code
                    plt.title(title)
                    plot_num += 1